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bispectrum estimators are expected to continue to provide the best constraints on the non-Gaussian parameters in 
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I. INTRODUCTION 



The standard inflationary paradigm predicts a flat Universe perturbed by nearly Gaussian and scale invariant pri- 
mordial perturbations. These predictions have been verified to a high degree of accuracy by Cosmic Microwave Back- 
ground (CMB) and Large-Scale Structure (LSS) measurements, such as those provided by the Wilkinson Microwave 
Anisotropy Probe (WMAP; [Komatsu et a/.| [2009)1, the 2dF Gal axy Redshift Survey (2dFGRS; [Percival et al.\^M2\ 



and the Sloan Digital Sky Survey (SDSS; Tegmark et al. 2004). Despite this success, it has proved to be difficult to 



discriminate between the vast array of inflationary scenarios that have been proposed by high-energy theoretical inves- 
tigations, or even to rule-out alternatives to inflation.. Since most of the present constraints on the Lagrangian of the 
inflaton field have been obtained from measurements of the two-point function, or power spectrum, of the primordial 
fluctuations, a natural step is to extend the available information is to look at non-Gaussian signatures in higher order 
correlators. 

The lowest order additional correlator to take into account is the three-point function or its counterpart in Fourier 
space, the bispectrum. Every model of inflation is characterized by specific predictions for the bispectrum of the 
primordial perturbations in the gravitational potential <l)(k). The bispectrum B(^(ki,k2,k^) of these perturbations is 
defined as 

<<D(ki)(D(k2)<D(k3)) = {Infdni^n^) B^ik.MM) , d-D 

where we have introduced the notation k,y s ki +k2 so that the Dirac delta function here is 60(^123) = ^oi^i +^2 + ^3)- 
Together with the assumption of statistical homogeneity and isotropy for the primordial perturbations, this implies that 
the bispectrum is a function of the triplet defined by the magnitude of the wavenumbers ki , k2 and k^ forming a closed 
triangular configuration. The current constraints that we are able to derive on the bispectrum B(s,(k\,k2,k-i) provide 
additional information about the early Universe; the possible detection of a non-vanishing primordial bispectrum 
in future observations would represent a major discovery, especially as it is predicted to be negligible by standard 
inflation. 

The cosmological observable most directly related to the initial curvature bispectrum is given by the bispectrum 
of the CMB temperature fluctuations, which provide a map of the density perturbations at the time of decoupling, 
the earliest information we have about the Universe. Current measurements of individual triangular configurations of 
the CMB bispectrum are, however, consistent with zero. Studies of the primordial bispectrum, therefore, are usually 
characterized by constraints on a single amplitude parameter, denoted by /nl, once a specific model for Bo is assumed. 
Since most models predict a curvature bispectrum obeying the hierarchical scaling Bi[,(k,k,k) ~ P^{k), with /'o(^) 
being the curvature power spectrum, the non-Gaussian parameter roughly quantifies the ratio /nl ~ B,p(k, k, k)/P^{k), 
defining the "strength" of the primordial non-Gaussian signal. In addition, we can write 

B^(kuk2,k,) = fyLF(kuk2,ki) , (1.2) 

where F{k\,k2,ki) encodes the functional dependence of the primordial bispectrum on the specific triangle config- 
urations. For brevity, the characteristic shape-dependence of a given bispectrum is often referred to simply as the 



bispectrum shape (a precise definition of the bispectrum shape function will be given in section II. A 1. Inflationary pre- 
dictions for both the amplitude /nl and the shape of Bo that are strongly model-dependent. Notice that the subscript 
"NL" stands for "nonlinear", since a common phenomenological model for the non-Gaussianity of the initial condi- 
tions can be written as a simple nonlinear transformation of a Gaussian field. Generically, of course, non-Gaussianity 
is associated with nonlinearities, such as nontrivial dynamics during inflation, resonant behaviour at the end of infla- 
tion ('preheating'), or nonlinear post-inflationary evolution. At the very least, future CMB and LSS observations are 
expected to be able to eventually detect the small last contribution. 

Perturbations in the CMB provide a particularly convenient test of the primordial density field because CMB temper- 
ature and polarization anisotropics are small enough to be studied in the linear regime of cosmological perturbations. 
Once the effects of foregrounds are properly taken into account, a non-vanishing CMB bispectrum at large scales 
would be a direct consequence of a non-vanishing primordial bispectrum. As we will see, while other CMB probes 
of primordial non-Gaussianity are available, such as tests of the topological properties of the temperature map based 
on Minkowski Functionals or measurements of the CMB trispectrum, the estimator for the non-Gaussian parameter 
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/nl has been shown to be optimal. We will focus mostly on this bispectrum estimator in the section of this review 
dedicated to the CMB. 

In the standard cosmological model, the large-scale structure of the Universe, that is, the distribution of matter 
and galaxies on large scales, is the result of the nonlinear evolution due to gravitational instability of the same initial 
density perturbations responsible for the CMB anisotropics. This is, perhaps, the most important prediction of the 
inflationary framework which provides a common origin for the CMB and large-scale structure perturbations as the 
result of tiny quantum fluctuations stretched over cosmological scales during a phase of accelerated expansion. The 
large-scale structure we observe at low redshift, however, is characterized by large voids and small regions with very 
large matter density, and it is therefore a much less direct probe of the initial conditions. The distribution of matter 
becomes a highly non-Gaussian field precisely as a result of the nonlinear growth of structures, even for Gaussian 
initial conditions. This non-Gaussianity is expressed, in particular, by a non-vanishing matter bispectrum at any 
measurable scale, including the largest scales probed by current or future redshift surveys. In this context, the effect 
of primordial non-Gaussianity, i.e. of an initial component in the curvature bispectrum, will constitute a correction 
to the galaxy bispectrum. It follows that the possibility of constraining or detecting this initial component is strictly 
related to our ability to distinguish it from other, primary sources of non-Gaussianity, that is the nonlinear gravitational 
evolution, and, in the case of galaxy surveys, nonlinear bias. 

The study of non-Gaussian initial conditions for large-scale structure has a relatively long history, with important 
contributions going back to the mid eighties. The standard picture that has been developed over the years, assumed 
that, at large scales, the effect of primordial non-Gaussianity on the galaxy distribution is simply given in terms of 
an additional component to the galaxy bispectrum. This is obtained, in perturbation theory, as the linearly evolved 
and linearly biased initial matter bispectrum, related to the curvature bispectrum Bg,{ki,k2,k^) by the Poisson equa- 
tion. Such component becomes subdominant as the gravity-induced non-Gaussian contribution grows in time. In this 
framework, as one can expect, high-redshift and large-volume galaxy surveys would constitute the best probes of the 



initial conditions. It has been shown, in fact, that proposed and planned redshift surveys, such as Euclid (Refregier 



\et al. \ |2010 [), should be able to provide constraints on the primordial non-Gaussian parameters comparable, if not 
better, than those expected from CMB missions such as Planck. What is more important, in the event of a detection 
by Planck, is that confirmation by large-scale structure observations will be required. 

Recent results from N-body simulations with non-Gaussian initial conditions, however, have revealed a more com- 
plex picture. The effect of primordial non-Gaussianity at large scales is not limited to an additional contribution to 
the galaxy bispectrum, but it quite dramatically affects the galaxy bias relation itself, that is, the relation between the 
matter and galaxy distributions. A surprising consequence is that it induces a large correction even for the galaxy 
power spectrum. Such an effect has attracted considerable recent attention and, remarkably, have placed constraints 
on the non-Gaussian parameter from current LSS data-sets which already appear to marginally improve on CMB lim- 
its. However, from a theoretical point of view, a proper understanding of the phenomenon remains to be properly 
developed and, for example, reliable predictions for the galaxy bispectrum are not yet available. Most importantly, 
as for general cosmological parameter estimation, a complete likelihood analysis aimed at constraining, or detecting, 
primordial non-Gaussianity in large-volume redshift surveys should involve joint measurements of the galaxy power 
spectrum and bispectrum, as well as possibly higher-order correlation functions. While we are still far from a proper 
assessment of what such analysis would be able to achieve, current results in this direction are very encouraging. 

This review is divided in four parts. In section |ll] we will first discuss initial conditions as defined in terms of 
the primordial curvature bispectrum and its phenomenology. We will then review the observational consequences 
of primordial non-Gaussianity on the CMB bispectrum, section III and on the large-scale structure bispectrum as 
measured in redshift surveys, section IV In both cases we will discuss theoretical models for the observed bispectra 
and technical problems related to the estimation of the non-Gaussian parameters, with the differences that naturally 
characterize such distinct observables. We also give an example of joint analysis using both CMB and large-scale 
structure when we consider the possibility of constraining a strongly scale-dependent non-Gaussian parameter /nl(^), 
emerging in some recently proposed inflationary models. 



k. 



FIG. 1 Triangle types contributing to the bispectrum corresponding to 'squeezed' or local configurations with ^3 ki, k2 (left), 
equilateral configurations with kj, x k[ x ki (centre) and flattened configurations with k^, x ki + ^2 (right). 



II. INITIAL CONDITIONS AND THE PRIMORDIAL BISPECTRUM 

A. The primordial bispectrum and shape function 

The starting point for this discussion is the primordial gravitational potential perturbation <l)(x, f) which was seeded 
by quantum fluctuations during inflation or by some other mechanism in the very early universe (f <K fdec)- When 
characterizing the fluctuations <1) we usually work in Fourier space with the (flat space) transform defined through 

(D(x, 0= r ;^e-'"''0(k,f). (II.3) 

The primordial power spectrum P(s,(k) of these potential fluctuations is found using an ensemble average, 

<0(k)<l)*(k')) = (27r)35o(k - k')Po« , (11.4) 

where we have assumed that physical processes creating the fluctuations are statistically isotropic so that only the 
dependence on the wavenumber remains, k - |k|. Recall that for nearly scale-invariant perturbations, the fluctuation 
variance on the horizon scale k ^ H is, almost constant A^^^ » k^ P,ii{k)l2n^ x const., implying P,s,{k) ~ k^^. 
The primordial bispectrum Bij,{k k2, k-i) is found from the Fourier transform of the three-point correlator as 

<(D(ki )0(k2)(l)(k3)) = {2nf 5o(ki23) B<i>(^i M,^). (II.5) 

Here, the delta function enforces the triangle condition, that is, the constraint that the wavevectors in Fourier space 
must close to form a triangle, ki -H ki -H ks = 0. Examples of such triangles are shown in fig.[T[ illustrating the basic 
squeezed, equilateral and flattened triangles to which we will refer later. Note that a specific triangle can be completely 
described by the three lengths of its sides and so, in the isotropic case, we are able to describe the bispectrum using 
only the wavenumbers ki,k2, k^. The triangle condition restricts the allowed wavenumber configurations iki,k2, kj,) to 
the interior of the tetrahedron illustrated in fig. [2] 

The most studied primordial bispectrum is the local model in which contributions from 'squeezed' triangles are 
dominant, that is, with e.g. ^3 <s: ki, k2 (as illustrated in the left of fig. [T]l. This is well-motivated physically as it 
encompasses 'superhorizon' effects during inflation when a large scale mode k^ (say) which has exited the Hubble 
radius exerts a nonlinear influence on the subsequent evolution of smaller scale modes ki, k2- Although this effect 
is small in single field slow-roll inflation, it can be much larger for multifield models. In a weakly coupled regime, 
the potential can be split into two components, the linear term Ol, representing a Gaussian field, giving the usual 
perturbation results plus a small local non-Gaussian term Onl ( |Salopek & Bond 1990[ l, 



0(x) = <1)l(x) -i- <Dnl(x) 

- <Dl(x) + /nl[<1>2(x)-<(I)2(x))], (II.6) 
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FIG. 2 Tetrahedral domain for allowed wavenumber configurations fci , ^3 contributing to the primordial bispectrum B(ki , k2, k^)). 
A regular tetrahedron is shown satisfying ki + k2 + kj < 2fcmax = ^K. 



where /nl is called the nonlinearity parameter. In Fourier space, the nonlinear term is then given by the convolution 

J3i. 



^nl(^) = /nl 



r 7^1>L(k + k')3>L(k')-(2;r)^5fl(k)<02) 



(11.7) 



From this we can infer, using ( |II.4| l, that the only non-vanishing contributions to the bispectrum ( |II.5[ ) take the form 

<0(ki)(l)(k2)<l>(k3)> = 2(2;r)35o(ki23) P^ih) PM . (II.8) 
Accounting for permutations, the local bispectrum becomes 



B^ihMM) = ifn^PMP^ih) + PMPM + PMPM] 



2/p 



NL 



A2 



3 

^1^2 j 



(11.9) 



Although this is a rather pathological function which diverges along the edges of the tetrahedron {i.e. when any 
ki 0), we can infer from it some basic properties of the bispectrum for any model which is nearly scale-invariant. 
For example, we can observe that the bispectrum at equal ki has the characteristic scaling. 



B^(k, k, k) = 2/nlA2 . 



(11.10) 



If we remove this overall k~^ scaling by multiplying (II.9i by the factor_(A:iA:2A:3)^, we note that on transverse slices 
through the tetrahedron defined by ^ s {kx + k2 + k^)/2 - const, (see fig.pj) the bispectrum only depends on the ratios 
of the wavenumbers, say ^2/^1 and kj,/ki. Indeed, it can prove convenient to characterize the bispectrum in terms of 
the following transverse parameters ( Fergusson & Shellard 2007 [ Rigopoulos et al. 2006a| l 



\(ki +k2 + h). 



a -(k2- h)lk , 



P^{k-ki)lk, 



(11.11) 



with the domains k < fcmax, < j6 < 1 and -{ I - P) < a < 1 -j6. The volume element on the regular tetrahedron of 
allowed wavenumbers then becomes dkidkidki, - k^dkdadp. 
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FIG. 3 Shape functions for the scale-invariant equilateral (left) and local (right) models, S(ki,k2,k2) = S (a, p) on transverse slices 
with 2^ = ^1+^:2 + ^3= const. 



These considerations lead naturally to the definition of the primordial shape function ( [Babich et al. 2004| l 



1 



S{kx,k2,h) = -^{kxk2hfB^{kx,k2,h) , 



(11.12) 

where is a normalization factor which is often chosen such that S is unity for the equal A;, case, that is, S (k, k, k) = I 
(we shall discuss alternatives to this rather arbitrary convention later). For example, the canonical 'local' model (II. 9 1 
has the shape 

k^ 



S'°^'\ki,k2,k,)^'^- 



k^ k^ \ 
1 _l_ 2 _|_ 3 

^2^3 kiki kxk2^ 



(11.13) 



Thus it is usual to describe the primordial bispectrum in terms of an overall amplitude /nl and a transverse two- 
dimensional shape S{ki,k2,k^) - S(&, yS), which incorporates any distinctive momentum dependence. Of course, if 
there is a non-trivial scale dependence, then the full three-dimensional dependence of S{kx,k2,k-}) on the fc, must be 
retained. 

There are other physically well-motivated shapes in the literature which have also been extensively studied. The 
simplest shape is the constant model 

S'°"'\kx,k2,k3)^l, (11.14) 

which, like the local model, has a large-angle analytic solution for the CMB bispectrum ( [Fergusson & Shellard[|2009| . 
The local model tends to be the benchmark against which all other models are compared and normalized, but for 
practical purposes the constant model is much more useful, given its regularity at both late and early times. The 
equilateral shape is another important case with (jBabich et al. 2004| l 

(kx +k2- k3)(k2 + k3- kxXki +ki - ^2) 



(11.15) 



^1^2^3 

While not derived directly from a physical model, it has been chosen phenomenologically as a separable ansatz for 
higher derivative models ( Creminelli 2003) and DBI inflation ( [Alishahiha et a/., ^2004) . The equilateral shape is 
contrasted with the local rnodel in fig.jsj 

Another important early result was the primordial bispectrum shape for single-field slow roll inflation derived by 



Acquaviva et a/.| ( |2003| l; |Maldacenaj2003) 



''(^1,^2,^3) ^ (3e- 277) 



^2^3 kxk^, 



kxk2 



+ e 



^ (6e-27?)5 



local 



{kx,k2,h) + 



equil 



{kxk\ + 5 perm.) + 4 
(^1,^2,^3), 



kik2ki 



(n.i6) 
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where e, 77 are the usual slow roll parameters. In the second line, we have noted that this shape can be accurately 



represented as the superposition of local and equilateral shapes. The coefficients in (II. I61, which include the scalar 
spectral index n — \ = —be + 2ri ^ -0.05, confirm that /nl ^ 1 and so standard single slow roll inflation cannot 
produce an observationally significant signal. Nevertheless, it is interesting to determine which shape is dominant in 
( II.16| l and to what extent other primordial shapes are independent from one another. 

Whether two different primordial shapes can be distinguished observationally can be determined from the correla- 
tion between the corresponding two CMB bispectra weighted for the anticipated signal to noise, as in the estimator 



(see next section) and the Fisher matrix analysis (see section III.H 1. However, direct calculations of the CMB bispec- 
trum can be very computationally demanding. A much simpler approach is to determine the independence of the two 



shape functions S and S ' from the correlation integral ( Fergusson & Shellard 2009 see also Babich et al. ( 2004 1) 



where we choose the weight function to be 



') = I S(ki,k2,k-i)S'{ki,k2,ki,) oj^iki , ^2, k-i)d^k , 

J'Vt 



(jj^(ki,k2,k2) 



1 



k\ + k2 + k^ 



reflecting the primary scaling of the CMB correlator The shape correlator is then defined by 



C(S,S') = 



F{S,S') 



^/F{S,S)F{S',S') 



(11.17) 



(11.18) 



(11.19) 



Here, the integral is over the tetrahedral region shown in fig. |2] taken out to a maximum wavenumber k < k^^^ 
corresponding to the experimental range / < f^nax for which forecasts are sought (with /'^ax ~ ^max)- The weight 
function cj^{k] , k2, k^) appropriate for mimicking the large-scale structure bispectrum estimator (see section IV.C.2i, 
would be different with varying scaling lawas introduced by the transfer functions for wavenumbers k above and below 
A;eq, the inverse comoving horizon at equal matter-radiation. Nevertheless, the l/k weight given in ( II.18| l provides a 
compromise between these scalings, and so shape correlation results should offer a useful first approximation. 

Below we will survey primordial models in the literature, showing how close the shape correlator comes to a full 
Fisher matrix analysis. However, here we note that the local shape (11.13 1 and the equilateral shape (III. 53 1 have 
only a modest 46% correlation. For the natural values of the slow roll parameters e x jj we find the somewhat 
surprising result that S^'^^'^ is 99.7% correlated with (and it cannot be easily tuned otherwise because 3e ~ rj is 
not consistent with deviations from scale-invariance favored observationally « - 1 < 0). Such strong correspondences 
are important in defining families of related primordial shapes, thus reducing the number of different cases for which 
separate observational constraints must be sought. 



B. General primordial bispectra and separable mode expansions 



The three shape functions ( 11.13 1, ( II.14 i and ( III. 53 1 quoted above share the important property of separability, that 
is, they can be written in the form 



S{kx,k2,h) = Xiki) Y{k2)Z(k'i) + 5 perms. 



(11.20) 



or as the sum of just a few such terms. As we shall see, if a shape S is separable, then the computational cost of 
evaluating the corresponding CMB bispectrum Be^t2ei is dramatically reduced. In fact, without this property, the task 
of estimating whether a non-separable bispectrum is consistent with observation appears to be intractable (for large 
^max)- Of course, the number of models which can be expressed directly in the form (11.20 1 is very limited, despite the 



usefulness of approximate ansatze such as the equilateral shape (III. 53 1. Indeed, approximating non-separable shapes 
by educated guesses for for the separable functions X, Y, Z is neither systematic nor computationally efficient (because 
arbitrary non-scaling functions create numerical difficulties, as we shall explain later). 



FIG. 4 The one-dimensional tetrahedral polynomials q„(k) on the domain pi.23^ , rescaled to the unit interval for n = 0-5. Also 
plotted are the shifted Legendre polynomials P„C^^ - 1) (dashed lines) which share qualitative features such as n nodal points. 



Instead, we shall present a separable mode expansion approach for efficient calculations with any non-separable 
bispectrum, as described in detail in Fergusson et al. ( 2009 1 (and originally proposed in Fergusson & Shellard 2QQ7\ . 
Our aim will be to express any shape function as an expansion in mode functions 



(11.21) 



where, here, for convenience we have represented the symmetrized products of the separable basis functions qpik) as 

Q„{ki ,k2,h)^\ [qp{x)qr{y)qs{z) + 5 perms] = q\p q,- q^} , (11.22) 

with a one-to-one mapping ordering the products as n «-» {prs}. The important point is that the qp(k) must be an inde- 
pendent set of well-behaved basis functions which can be used to construct complete and orthogonal three-dimensional 
eigenfunctions on the tetrahedral region defined by (see fig.|2]) 

ki , k2, k^ < kmax , ki < k2 + ky, for k\ > k2, kj,, + 2 perms . (11.23) 

The introduction of the cut-off at fc^ax is motivated both by separability and the correspondence with the observational 
domain I < £max- In the shape correlator (11.19 1, we have already seen what is essentially an inner product between 
two shapes on this tetrahedral region, which we can define for two functions /, g as 



r f(ki,k2, k3)g(kuk2, ki) co(ki,k2, h)d^r ■. 

J'Vt 



(11.24) 



with weight function w. 

Satisfactory convergence for known bispectra can be found by using simple polynomials qp(k) in the expansion 
( |II.21| i, that is, using analogues of Legendre polynomials on the domain (11.23 1. With unit weight, the polynomials 



satisfying {qp{k\), qr{k\)) — dp,- can be found by generating functions with the first three given by (Fergusson et al. 
[20091 ) 



^o(x)=V2, ^i(x) = 5.79(-i^+x), ^2(x) = 23.3( 



54 

215 43 



4«X + x2) 



(11.25) 



The first few polynomials qp{k) are plotted in fig.|4] where they are contrasted with Legendre polynomials. 

The three-dimensional separable basis functions Q„ in (11.22 1 reflect the six symmetries of the bispectrum through 
the permuted sum of the product terms. They could have been constructed directly from simpler polynomials, such as 
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FIG. 5 Orthonormal eigenmode decomposition coefficients l |II.26[ l for the equilateral and DBI models (upper panel) and shape 
correlations pi.l9[ ) of the original bispectrum against the partial sum up to a given mode n (lower panel). The correlation plot 
includes both primordial and late-time CMB bispectra for the equilateral and DBI models, as well as the late-time CMB bispectrum 
from cosmic strings (refer to section [Hill. In all cases, we find that we need at most 15 three-dimensional modes to obtain a 
correlation greater than 98% (primordial convergence without the acoustic peaks requires only 6 modes). 



1 , ki +k2 + k2, k^ + k^+k^, however, the qp polynomials have two distinct advantages. First, the ^^'s confer partial 
orthogonality on the <3„ and, secondly, these remain well-behaved when convolved with transfer functions. 

In order to rapidly decompose an arbitrary shape function S into the coefficients aJJ «-> a^^j, it is more convenient 
to work in a non-separable orthonormal basis ((!??„, %„) - 6„,„. These can be derived directly from the Q„ through 



Gram-Schmidt orthogonalization, so that "R„ = Tj'p=q ^mpQp with A^p a lower triangular matrix (see Fergusson et al. 
[2009| ). Thus we can find the unique shape function decomposition 

N N N 

S(kuk2,k3) = ^<K(fci,^2,^3) = 2''«^«(^i'^2,^3), With a« = {S , K) and = 2(^^)„^a« . (11.26) 

n n p 

In the orthonormal frame, Parseval's theorem ensures that the autocorrelator is simply {S, S) - 2„ a'f'. Hence, 
with a simple and efficient prescription we can construct separable and complete basis functions on the tetrahedral 



domain ( 11.23 1 providing rapidly convergent expansions for any well-behaved shape function S . These eigenmode 
expansions will prove to be of great utility in subsequent sections. Examples of this bispectral decomposition and its 
rapid convergence for the equilateral and DBI models are shown in fig. [5] 
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C. Families of primordial models and their correlations 



We will now briefly survey the main categories of primordial models in the literature and their relative independence, 
closely following the discussion in jFergusson & Shellard ( 2009| l. 



1 . The constant model 



The constant model ( II. 14 1 is the simplest possible primordial shape with triangles of every configuration contribut- 
ing equally to the bispectrum ^2, ^3); it is the equipartition model. The constant model was motivated initially 
by its simplicity (|Fergusson & Shellard 2009 1 leading to an analytic solution for the large-angle CMB bispectrum, as 
well as due to its close correlation with equilateral models. However, the shape does have more explicit physical moti- 



vation in at least one context (Chen & Wang 2009j), during multifield inflation for a slowly turning trajectory (denoted 
quasi-single field inflation). For multifield inflation, it is well known that the conversion of isocurvature fluctuations 
into curvature fluctuations during 'corner-turning' can source significant non-Gaussianity (see e.g. Rigopoulos et al. 
2006a[[\irmzzi & Wands|[2006^ . In the quasi-single field case with mass m ~ H isocurvature modes, a detailed inves- 



tigation of the ongoing conversion into the curvature mode demonstrated that novel shapes could be generated ( Chen & 



Wang 2009 1, amongst them shapes which were very nearly constant. Generically, these model-dependent shapes be- 
longed to a one-parameter family which interpolated non-trivially between equilateral ( |IlI.53| l and local ( |n.l3| l shapes 
(see also Renaux-Petel 2009 Senatore et al. 2009 1. This is an important caveat for the present discussion, because 



non-Gaussian searches could uncover shapes intermediate between the categories we will discuss below. 



2. Equilateral triangles - centre-weighted models 

Bispectra dominated by contributions from nearly equilateral triangle configurations, k\ ^ k2 ~ can be fairly 
easily characterized analytically and are the most amenable to CMB searches. However, equilateral non-Gaussianity 
requires the amplification of nonlinear effects around the time modes exit the horizon, which is not possible in a slow- 
roll single field inflation. Instead, the kinetic terms in the effective action must be modified as in the Dirac-Born-Infeld 
(DBI) model ( Alishahiha et al. 2004)l or by explicitly adding higher derivative terms, such as in K-inflation (see, for 
example, Chen et al. 2007 | l. The resulting corrections modify the sound speed and inflation is able to take place 



in steep potentials. For DBI inflation, this leads to non-Gaussianity being produced with a shape function of the form 
( fAlishahihaTFari [20041 [Creinineni, ,2003 J 



S(ki,k2,ki) - 



1 



k\k2kT,{k\ + k2+ k^)^ 



(11.27) 



Another example of a model with non-standard kinetic terms is ghost inflation ( Arkani-Hamed et al. 2004 1 with a 
derivatively-coupled field driving inflation and a trilinear term in the Lagrangian creating a non-zero equilateral-type 
shape S tending towards constant. 

General non-Gaussian shapes arising from modifications to single field inflation have been extensively reviewed in 
Chen et aL\ ( |2007[ ). Using a Lagrangian that was an arbitrary function of the field and its first derivative, they were 
able to identify six distinct shapes describing the possible non-Gaussian contributions. Half of these had negligible 





DBI 


Ghost 


Single 




C(5,5') 


C(S, B') 


C(5,5') 


C(B,B') 


C(5,5') 


C(B,B') 


Equilateral 


0.99 


0.99 


0.98 


0.98 


0.95 


0.96 


DBI 






0.94 


0.95 


0.98 


0.99 


Ghost 










0.86 


0.89 



TABLE I Shape correlations ( |II.19| l and CMB correlations ( |III.61^ between the equilateral family of primordial models. 
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Local 


Warm 


Flat 


Feature 




C(5,5') 


C(S, B') 


C(5,5') 


C{B,B') 


C(5,5') 


C(B,B') 


C(5,5') 


C(B,B') 


Equilateral 


0.46 


0.51 


0.44 


0.42 


0.30 


0.39 


-0.36 


-0.43 


Local 






0.30 


0.52 


0.62 


0.79 


-0.41 


-0.39 


Warm 










0.01 


0.21 


-0.05 


-0.27 


Flat 














-0.44 


-0.32 



TABLE II Shape correlations jll. 19[ > and CMB correlations l |III.61[ l for 5 distinct families of primordial non-Gaussian models. 



amplitude being of the order of slow-roll parameters (with two already given in ( 11. 16 1). Of the remaining three shapes 
( Chen etal. 2007 see also Seery & Lidsey 2005) , one was believed to be subdominant, the s econd recovered the DBl 
shape ( 11.27 1, leaving a third distinct single field shape which is the inverse of the local shape ( 11. 1 3 i, 5 ^'"^''^ oc S ' . 
Finally, we recall the original equilateral shape (111.53 i, noting that it was introduced not because of a fundamental 
physical motivation, but as a separable approximation to the DBl shape ( 11.27 1 ( Babich et al. |2004 1. 

Despite the apparent visual differences between these shapes (see |Fergusson & Shellard| |2009 , particularly near 
the edges of the tetrahedral domain, the shape correlator (11.19|l reveals at least a 95% or greater correlation of the 



DBl, ghost and single shapes to the equilateral shape (111.53 i (consistent with results in Babich et al. 2004 Smith 



& Zaldarriaga 2006 1 . Comparative results between the shape correlator are given in Table |I| (together with the 



corresponding CMB correlation results brought forward and showing the efficacy of these estiomates). These particular 
centre-weighted shapes must be regarded as a single class which would-be extremely differentiate observationally, 
without a bispectrum detection of very high significance. 

Finally, we comment on the 'orthogonal' shape 5™*og proposed inll 



Senatore et al. 



( |2009|, together wi th 5'=''"", for 
characterizing single field inflation models with an approximate shift symmetry (see also Chen et al. 2007 1. This shape 
is approxim ately ^Q'^og oc ^^qmi _ 2/3, which means it is very similar to an earlier study of flattened shapes (Meerburg 
et al. 2009 ) which proposed an 'enfolded' shape with 5''="f°i'* oc ^^quii _ p^om the eigenmode decomposition (11.26 1 



of the equilateral model shown in fig.|5] it is clear how the degree of correlation can be altered by subtracting out the 
important constant term. 



3. Squeezed triangles - corner-weighted models 



The local shape covers a wide range of models where the non-Gaussianity is produced by local interactions. These 
models have their peak signal in "squeezed" states where one ki is much smaller than the other two due to non- 
Gaussian ity typ ically being produced on superh orizon scales. We hav e already observed t hat single-field slow-roll 



IP - - - , , 

inflation (|ll.l6l is dominated by the local shape (Bartolo et al. 2004a I, though f^^- is tiny (Acquaviva et al. 

[Baitolo eFoT, ,2004a, Maldacenaj 



E003 



). The production of non-Gaussianity during multiple held mnation 



2003 



(|Bernardeau & Uzan||2003||2002[|Lyfli & Rodriguez] [2005] [Rigbpoulos et a/.||2006a|b| [20071 [Seery & Lidsey[[2005 



Vernizzi & Wands 2006 1 shows much greater promise through conversion of isocurvature into adiabatic perturbations 



(see, for example, recent work in Byrnes et al.\ [2008[ [Chen & Wang , 2009] [Naruko & Sasaki[ [2009, Renaux-Petel] 
and references therein). The magnitude of the non-Gaussianity generated is normally around ~ 0(1), which 
is at the limit for Planck detection, but models can be tuned to create larger signals. Significant f^^' can be produced 



2009 



in curvaton models with f^^- ^ (9(100) ( ^Bartolo et al. 



2004b 



Linde & Mukhanov, 



2006;,Lyfli et al. 



f^- can also be generated at the end of inflation from massless preheating or other reheating mechanisms ([Chambers 



20031. Large 



& Rajantie[[2008[|Enqvist et a/.[[2005a[ b') 

We note that local non-Gaussianity can also be created in more exotic scenarios. Models based on non-local field 
theory, such as p-adic inflation, can have inflation in very steep potentials. Like single field slow-roll inflation, the 



predicted 'non-local' shape function is a combination of a dominant local shape ( [11.13| l and an equilateral shape (111.53 1 



(see, for example,|Barnaby & Cline. 2008 ; Halliwell, 1993 ; Seery etal. 2008 ; Zaldarriaga 2004). The e kpyrot ic model 



can also generate significant f^^- ([Buchbinder et al. 2008 ; Creminelli & Senatore 2007 ; Koyama et a/.[[2007[[Lehners[ 



& Steinhardt 2008a|b[ ). Here the density perturbations are generated by a scalar field rolling in a negative exponential 



potential, so non-linear interactions are important with f^-^- x (9(100). 
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In using the shape correlator for the local model, we must introduce a small-wavenumber cut-off, taken to be a 
^min = 2/tau(), otherwise the shape correlator C(5'°'^''', 5'°'^^') becomes infinite. This logarithmic divergence does not 
afflict the CMB bispectrum because we do not consider contributions below the quadrupole Z = 2 (a threshold 
which is approximated by the primordial cut-off). The local shape is modestly correlated at the 40-55% level with 
the equilateral shapes, mainly through the constant term in the expansion ( |1I.26[ ). As can be seen in table |ll] this 
somewhat underestimates the CMB correlator Nevertheless, a NG signal of only modest significance should be able 
to distinguish between these independent models. 

Finally, warm inflation scenarios, i.e. models in which dissipative effects play a dynamical role, are also predicted 
to produce significant non-Gaussianity ( [Gangui et al. 1994 Gupta et al. 2002 1. Contributions are again dominated 
by squeezed configurations but with a different more complex shape possessing a sign flip as the comer is approached 
(see fig.|6]l. This makes the warm 5™™ and local 5'i°<=^' shapes essentially orthogonal with only a 33% correlation. 
Again, in using the shape correlator, we need to introduce the same phenomenological cut-off k^in as for the local 
model, but we also note the more serious concern which is the apparent breakdown of the approximations used to 
calculate the warm inflation shape near the corners and edges. 



4. Flattened triangles - edge-weighted models 



It is possible to consider inflationary vacuum states which are more general than the Bunch-Davies vacuum, such as 



an excited Gaussian (and Hadamard) state (Holman & ToUey 2008 see also discussions in Chen et al. 2007 Meerburg 
et g/. 1 12009) . Observations of non-Gaussianity in this case might provide insight into trans-Planckian physics 



The 



proposed shape for the bispectrum is 



+ 2 perms \ + 



2{k\ + kl + kl)s 



(ki +k2- k^nki +h- +h- k2)^ 



(11.28) 



The bispectrum contribution from early times is dominated by flattened triangles, with e.g. k-i ^ k\ +k2, and for a small 
sound speed Cj <sc 1 can be large. Unfortunately, as the divergent analytic approximation breaks down at the boundary 
of the allowed tetrahedron, some form of cut-off must be imposed, as shown for the smoothed shape in fig. |6] where 
an edge truncation has been imposed together with a Gaussian filter The lack of compelling physical motivation and 
ill-defined asymptotics make predictions for this model uncertain. 
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5. Features - scale-dependent models 



There are also models in which the inflation potential has a feature, providing a break from scale-invariance. This 



can take the form of a either a step (Chen et al. 20061 or a small oscillation superimposed onto the potential (Bean 



et al. 2008|. Analytic forms for both these three point functions have been presented in Chen et al. (2008 1 with one 
approximation taking the form, 



S^'^'ihMM) sin^ 



k\ + k2 + ki 



+ P 



(11.29) 



where k* is the associated scale of the feature in question and f is a phase factor Results for the shape correlator for a 
particular feature model (with k* « £* It^ and (* - 50), are given in table |n] showing that it is essentially independent 
of all the other shapes. Clearly, scale dependent feature models form a distinct fifth category beyond equilateral, local, 
warm and flat shapes. 



III. COSMIC MICROWAVE BACKGROUND 
A. The CMB bispectrum 

In this section we will study the connection between the primordial bispectrum at the end of inflation and the 
observed bispectrum of CMB anisotropics Bf^e^c^. Our work will be primarily concerned with the analysis of the 
three-point function induced by a NG primordial gravitational potential C>(k) in the CMB temperature fluctuation 
field. Temperature anisotropics are represented using the a(,„ coefficients of a spherical harmonic decomposition of 
the cosmic microwave sky, 

AT T-, 7- 
—in) ^ 2^ai,Jtm(n) . 

Em 

Analogous expansions are performed for the E mode polarization field in order to produce polarization multipoles a^^. 
For simplicity and clarity, throughout most of this review we will focus on the temperature multipoles a^^^^, and omit 
the superscript T for convenience of notation. However we stress here that all of the considerations we make in the 
following can be readily applied to polarization multipoles and related bispectra. More discussion about this subject 
can be found in section llIl.H.2l 

The primordial potential d) is imprinted on the CMB multipoles «/,„ by a convolution with transfer functions Ai(k) 
representing the linear perturbation evolution, through the integral: 

flf,„ = An{-i)' r Afik) (D(k) Ff„,(k) . (III.30) 

The radiation transfer functions encode all the typical effects observed in the CMB power spectrum at linear order, 
that is, the Sachs-Wolfe effect. Integrated Sachs-Wolfe effect, doppler peaks and silk damping (see e.g. |Dodelson] 
2003 Ma & Bertschinger 1995). An equation identical to ( 111.30| l produces the E-mode polarization CMB multipoles 



starting from the primordial temperature fluctuation field, provided polarization transfer functions replace temperature 



transfer functions in the convolution above. It is sometimes useful to rewrite equation (III. 30 1 in position, rather than 
Fourier, space. In this case it is straightforward to show that ( |III.30| becomes 



drr^ae{r)^tm{r) , (III.31) 



where, starting from the primordial potential 0(x), we transform from Cartesian into polar coordinates x - (r, x), and 
defined: 

^fjr) = I £/Q,<l)(x)F,„(x) , (III.32) 



^ Ij dkk^Ae(k)jt(kr) . (III.33) 
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In this expression j( is the spherical Bessel function of order {. The CMB bispectrum is the three point correlator of 
the ac„,, so substituting we obtain 



Bmmiim - ko-iMn^acim.ahmz) (111.34) 

= (4;r)^(-0'>^'^^'^ r |^H^A,fc)A,fe)A,fe)x (111.35) 
J {Iny {Inf (2ny 

<0(ki)0(k2)<D(k3)) Ff„„,(ki) i-f.^.d^z) Y(,„,(h) (111.36) 

- J^J J X^t/x J J^it!!fe2fl?^3(fel^2fe3)'B4.(^I,^2,fe3)Af,(^l)Af,(^2)Af3(^3) (111.37) 

Xjc,(kix)je,(k2x)j{^(k2x) J c/Qj Ff,„„(x)F^j,„,(x)F^„„3(x) , (111.38) 

where in the last line we have integrated over the angular parts of the three k,, having inserted the exponential integral 
form for the delta function in the bispectrum definition ( 11.5 1. The last integral over the angular part of x is known as 
the Gaunt integral, which can be expressed in terms of Wigner-3 j symbols as (for more details on these functions and 
their properties (see e.g. |Komatsu] |2002| and references therein) 



'(2{i + l)(2(2 + 1)(2^3 + 1) I ti £2 i^M ti h (2 



4ji \ / \ mi m2 nij, 



(111.39) 



Given that most theories we shall consider are assumed to be isotropic, the m-dependence can be factorized out of the 
physically relevant part of the bispectrum ( |Luo||1994 l. It is then usual to work with the angle-averaged bispectrum, 

) . (111.40) 



/—i \ mi 

m: > 



m2 niT, 



or the even more convenient reduced bispectrum which removes the geometric factors associated with the Gaunt 
integral, 

gtihh ^Qlxhh L (11141) 

From the previous two formulae we also derive the following useful relations between the full, averaged and reduced 
bispectra: 



(2^1 + 1)(2^2 + 1)(2^3 + 1) / ^1 i2 u>rnm2m, _ i d h ^3 



The reduced bispectrum from ( 111.34| i then takes the much simpler form 



dkidk2dki (^1^2^:3)^ B,^(k\,k2,k-i) 

X Af , (^1 ) (k2) Ae, ih) je, (k, x) (k2x) je, ihx) . (111.43) 

This is the key equation in this section, since it exphcitly relates the primordial bispectrum, predicted by inflationary 
theories, to the reduced bispectrum observed in the cosmic microwave sky. This formula is entirely analogous to the 
well known relation linking the primordial curvature power spectrum P(s,{k) and the CMB angular power spectrum Q, 
that is. 



Cf = - j dkk^P^(k)Ajik) . (111.44) 

Finally, it is important to note that the Gaunt integral in ( 111.4 l| i encodes several constraints on the angle averaged 
bispectrum B(,e^(^ which are no longer transparent in the reduced bispectrum bi^t^e-^. These are: 
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1. The sum of the three multipoles ^, must be even. 

2. The ti's satisfy the triangle condition \£i - (j\ < 4 < 4 + ij- 

Analogous to the wavenumber constraint ( |11.23| l, the second condition is tells us that the only multipole configurations 
giving non-zero contributions to the bispectrum are those that form a closed triangle in harmonic (^-)space. For 
wavenumbers, the triangle condition is enforced through the x-integral over the three spherical Bessel functions jc(kix) 
which evaluates to zero if the A;, 's cannot form a triangle, whereas in multipole space it is enforced by the angular 



integration dQ.x over the spherical harmonics in (III. 39 1 



B. Separable primordial shapes and ClUIB bispectrum solutions 



In terms of the shape function (11.12 1, the reduced bispectrum (111.43 i can be rewritten as 

btihti = ^(^) \ x^dx J dkidk2dk3Sikuk2,ks)A{,iki)A{^ik2)A{j(ki)j{j(kix)je^ik2x)jej(k3x). (III.45) 

The expression above can be simplified and simple analytic solutions can sometimes be obtained for the very important 
class of separable shapes obeying the ansatz S - XYZ, as in ( 11.20 1. Substituting ( 11.20 1 into ( III.45 1, we find that 



I 



drr (r) (r) Zc, (r) + 5 perms 



(III.46) 



where we have defined the quantities: 



Ze(r) 



dkk^X(k) je(kr)Ae 



dkk^ Y(k) jiikr) At 



dkk^Z(k) ji(kr)Ae 



(III.47) 



Instead of the three-dimensional integral of ( III.45 1 we now have to deal with a much more tractable product of three 
one-dimensional integrals. Moreover, if we work at large angular scales in the Sachs-Wolfe approximation, the transfer 
functions become Ai{k) - \ji [(t,, - Tdec)k]. The presence of a product of spherical Bessel functions in the integrals 
above can lead in some cases to simple analytic solutions. 

Let us demonstrate this for the separable primordial shapes considered in sectionjll] The simplest possible shape, the 



constant model ( II. 14 1 with S{k\,k2,kj,) - 1 , has a large-angle analytic solution for the reduced bispectrum ( .Fergusson 
|&ShelEd| 120091 1, 



1 



HN(ltx H- l)(2^2 + 1)(2^3 + 1) 



1 



1 



ti+h + h + 'i £i+(2 + h 



(I <K 200) . 



(III.48) 



The large-angle solution ( |III.48| l is an important benchmark with which to compare the shape of late-time CMB 
bispectra from other models bf^(^Ci (note the /""^ scaling). The more general constant solution does not have an analytic 
solution because the transfer functions cannot be expressed in a simple form, but it can be evaluated numerically from 
the expression 



onst _ j 



x^dxl{^{x)l(^{x)li^{x) , where Ii(x) 



dk Ae(k) jt(kx) . 



(III.49) 



The numerical solution is shown in fig.|7] exhibiting the a regular pattern of acoustic peaks introduced by the oscillating 
transfer functions. 
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FIG. 7 The reduced CMB bispectrum for the constant model fcf™ff^ normalized relative to the large-angle constant solution i III.48 1. 
On the left, the bispectrum is plotted over the allowed tetrahedral region (see fig.[2]l using several density contours (light blue positive 
and magenta negative) out to f, < 2000 and, on the right, transverse triangular slices are shown at ^1+^2 + ^3= 2000. Note the 
coherent pattern of acoustic peaks induced by the transfer functions. From jFergusson et a/.H2009| l. 



For the local shape ( II.13| l, the Sachs-Wolfe approximation also yields a large-angle analytic solution 

1 1 1 



[^local 



2A^ 



277r2\A(A + 1X2(^2 + 1) ^2(^2 + 1X3(^3 + 1) ^3(^3 + lXi(A + 1) 



(III.50) 



where the divergences for the squeezed triangles (^1 <k k2,k-i...) in the primordial shape (11.13 1 are also reflected in 



b'f'^f'l . It is straightforward, in principle, to calculate the full bispectrum from the separable expressions arising from 
( |ffl3l l. 



(III.51) 



where the separated integrals analogous to (III.49 1 become 

aiix) J dkk^ A((k) itikx) , Pfix) J dkk^P^(k) A({k) j({kx) . (III.52) 
However, we note that these highly oscillatory integrals must be evaluated numerically with considerable care. 



For the equilateral shape (111.53 1 we first make its separability explicit by expanding the expression in the form, 

' k^ 



S(ki,k2,ki) 



2 \ Ik \ 
- — 1-2 perm. + 1- 5 perm. 



(III.53) 



While there is no simple large-angle analytic solution known for the equilateral model, it can be evaluated from the 
simplified expression 

^^"k ^ J ^^^^ [26(,d(Jc, + [a{^j3(^/3(j + 2 perm.) + [Pcjc^Se, + 5 perm.)] , (III.54) 



where cc/, j8, are given in (III. 52 1 and 7,, 5i are defined by (compare with the local case) 



7i{x) 



dkk^ p^(ky'^Ai(k) ji(kx) 



Si(x) 



dkk^P^ikfl^Aiik) ji(kx). 



(III.55) 
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C. Non-separable bispectra revisited 



Recall the mode expansion ( n.21|i of a general non-separable primordial shape. If we substitute this into the expres- 
sion for the reduced bispectrum (111.45 1, then the separability of the expansion leads to the same efficient calculation 
route discussed in the previous section through ( jFergusson et al. 2009 | l 



'Yfp'-' J^x^t/xj ^ ^ dhqp{k\)^c^(ki) ii^ikyx) J dk2qr(k2) Ae.iki) je.ikix) 



prs ^ 



dk-i qsih) Ae^iki) je^ikjx) 



+ 5 perm. 



(III.56) 



2 J ti e-, t 
X dxq^'^q^-q^ 



where the q^^ simply result from convolving the basis functions qp{k) with the transfer functions, 

= \ J dkqp(k)A((k)j((kx) . 



(111.57) 



The computationally costly 3D integrals have again reduced to a sum over products of ID integrals; we note that this 
economy arises because the triangle condition is enforced in ( 111.56 1 through the product of Bessel functions, resulting 
in a manifestly separable form in which we can interchange orders of integration. With this mode expansion, all non- 
separable theoretical CMB bispectra bi^e^e, become efficiently calculable provided there is a convergent expansion for 
the shape function. 

In the same way that we decomposed an arbitrary primordial shape S(ki,k2,kj) in section 11. B it is possible to 
construct analogous late-time separable basis functions Q„ and orthonormal modes with which to describe the 
CMB bispectrum fSf^f^ji (Fergusson et a/. , 2009; Fergusson & Shellard 2007). The tetrahedral domain defined by 
the triangle condition for multipole configurations {€1,^27^3} is essentially identical to that for wavenumbers (11.23 1, 
except that only even cases contribute 2 ^1 -H ^'2 + ^3 - 2n, n e N. However, the appropriate weight function now 
incorporates Wigner-3 j symbols arising from bispectrum products, 

\2 



1 

4^ 



(2^1 + 1)(2^2 + 1)(2^3 + 1) 



/i h h 




and 



~2 2 2~ 



(111.58) 



where in the second expression we have exploited the freedom to divide by a separable function V( - {2( + 1)'^^ and 
use a weight which makes the bispectrum functions more scale-invariant (eliminating an ^"'^^ factor - see below). The 
inner product between two functions /fif^f, and gt^t^c^ is altered from the primordial wavenumber integral (??) into a 
sum over multipoles on the tetrahedral domain, that is. 



(111.59) 



But for the change in the weight (which only affects configurations near the edges of the tetrahedron), the ID poly- 
nomials qp{€) and 3D separable product basis functions Q„{(i,(2, ^3) - ^{pQrls} (n <-» {prs}), as well as the resulting 
orthonormal modes are nearly identical to their primordial counterparts qpik), Q„(kuk2,kj,) and 'R„{ki,k2,k3) de- 
fined in section HLB] 

We can now expand an arbitrary QVIB bispectrum beiC^e, in both the separable and orthonormal mode expansions, 
which is achieved in the following form. 



(111.60) 
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where the variance term y]CcC(Cc reflects the signal-to-noise weighting expected in the CMB estimator (see sec- 
tion III. El. Again, the coefficients in the expansions are determined, first, from the orthonormal in ner pro ducts 
ff" = ■) and, secondly, the separable ajj are found with the transformation matrix analogous to ( 11.26 1. Ex- 
amples of the convergence of these mode expansions for equilateral, DBl and cosmic string CMB bispectra are given 
infig.|5] 



D. CMB bispectrum calculations and correlations 



Prior to the systematic mode expansion approach (111.56 1 being implemented, robust hierarchical schemes were de- 
veloped to calculate any non-separable CMB bispectru m ([111.45 1 directly (Fergusson & Shell ard , 200 7 2009). These 



use the transverse coordinate system {k, a, 0) given in ( 11. 11 1 and employ adaptive methods on a triangular grid to 



accurately determine the oscillatory 2D o-yS-integrations, with important efficiencies also coming from the flat sky ap- 
proximation, binning and interpolation schemes. Precision to greater than 1% across the full Planck domain f < 2000 
was established by direct comparison with analytic solutions such as (111.48 1 and (111.50 1. Examples of nonseparable 
(and separable) CMB bispectra found using these hierarchical coarse-graining methods are shown in figs [7] and [8] 
While the CMB bispectra be^f^j^ retain the qualitative features of the primordial shape functions S(k\,k2,k-i), they are 
overlaid with the oscillatory transfer functions which give rise to a coherent pattern of acoustic peaks. These direct 
bispectrum calculations revealed that typical primordial models could be described by eigenmode or other expansions 
using only a limited number of terms. 

Motivated by the form of the CMB estimator, we can define the following correlator to determine whether or not 
two competing theoretical bispectra can be distinguished by an ideal experiment. 



CcCr^Ce, 



(III.61) 



where the normalization is defined as 



Z 



d2 



D/2 



(111.62) 



The emergence of the inner product (111.59 1 in the expression (111.61 1 means that substitution of the mode expansions 
( |111.60[ ) for the theoretical bispectra reduces the correlator to 



(111.63) 



While the late time correlator (111.61 1 is the best measure of whether two CMB bispectra are truly independent, it can 



be demonstrated that for the majority of models the shape correlator ( 11. 19 1 introduced earlier is sufficient to determine 
independence. 

On the basis of the direct calculation of the bispectrum results and the CMB correlator, we can now quantitatively 
check the forecasting accuracy of the primordial shape correlator proposed previously (again closely following the 
discussion in Fergusson & Shellard (2009|l). 



1 . Nearly scale-invariant models 

For nearly scale-invariant models, the centre values for the bispectrum bm all have roughly the same profile but 
with different normalisations. As we see from (??), the oscillatory properties of the transfer functions for the CMB 
power spectrum, create a series of acoustic peaks for any combinations involving the following multipole values, 
I - 200,500,800,.... Of course, to observe the key differences between the scale invariant models we must study 
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the bispectmm in the plane orthogonal to the (I, I, Z)-direction, that is, the directions reflecting changes in the primor- 
dial shape functions. To plot the bispectrum (see figs |7] and [8]l, we consistently divide bm by the large-angle CMB 
bispectrum solution for the constant model (II.14i. This is analogous to multiplying the power spectrum Q's by 
1(1 H- 1), because it serves to remove the overall {^'^ scaling of the bispectrum, flattening while preserving the transverse 
momentum-dependence primordial shape and the effects of the oscillating transfer functions. 

The starting point is the constant model (11.14 1 which, despite its apparent simplicity, has a CMB bispectrum b'j°(^l 
revealing a non-trivial and coherent pattern or acoustic peaks that we have already noted (see fig. |7]). Given that 
the constant model has no momentum dependence, we stress that the resulting bispectrum is the three-dimensional 
analogue of the angular power spectrum {({+l)Cc for a scale -invariant model. The largest (primary) peak, for example, 
is located where all three fj - 220 (corresponding to the large blue region near the origin). We can interpret fig. |7] 
therefore, as the pure window function or beam effect of convolving any model with the radiation transfer functions 
A({k) while transforming from Fourier to harmonic space. 

The CMB bispectrum for the equilateral model is plotted in fig. |8] showing how the the centre-weighting from 
the primordial shape is well-preserved despite the convolution with the oscillating transfer functions. For the full 



CMB correlator (III.61 1, the DBI, ghost and single shapes are generally even more closely correlated with equilateral, 
presumably because distinctive features are 'washed out' by the transformation from Fourier to harmonic space. Com- 
parative results between the shape correlator and Fisher matrix analysis are given in Table |lj establishing that these 
models are highly correlated and diflicult to set apart observationally. 

The CMB bispectrum for the local model is also shown in fig. [8] demonstrating a marked contrast with equilateral 
which reflects their different primordial shapes shown in fig. [3] The dominance of the signal in the squeezed limit 
creates strong parallel ridges of acoustic peaks which connect up and emanate along the comer edges of the tetrahe- 



dron(see Bucher et al. 2009 , for further details). The 51% CMB correlation between the local and equilateral models 
is underestimated by the shape correlator at 41%, presumably because of effective smoothing due to the harmonic 
analysis. Reflecting their distinctive primordial properties, the CMB bispectra for the flat and warm models are poorly 
correlated with most of the other models, though the flat shape could be susceptible to confusion with the local CMB 
bispectrum with which it has a larger correlation (see table It is clear that the local, equilateral, warm and flat 
shapes form four distinguishable categories among the scale-invariant models. 



2. Scale-dependent models, cosmic strings and other late-time phenomena 



Models which have a non-trivial scaling, such as the feature models, can have starkly constrasting bispectra as 
illustrated in fig. ??. For example, instead of having the same pattern of acoustic peaks which characterise the scale- 
invariant models, the feature model can become entirely anticorrelated so that the primary peak has the opposite 



sign. Later, for this particular choice of k* in (11.29 1), for increasing / the phase of the oscillations becomes positively 



correlated by the second and third peaks. This can lead to small correlation with the other primoridial shapes, all below 
45% as shown in table[ll]for this k* and {nvdx- Clearly, these non-separable feature models form a distinct fifth category 
beyond the four scale-invariant shapes noted above and, of course, there are many possible model dependencies which 
can lead to further subdivision. 

By way of further illustration of the breadth of other possible non separable CMB bispectra, we present the late-time 
CMB bispectrum predicted analytically for cosmic strings (Regan & Shellard, ,2009) 



, stnng 
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where 4„„ = min(^i,^2,^3), t* = min(500, ^.m), C = min( 1/500, l/4,m) and 



{( < 2000) , 



(III.64) 



(III.65) 



Here, A ~ (^nGfj)^ is a model dependent amplitude with G/u - lu/mh measuring the string tension /u relative to the 



Planck scale. The cutoffs around i x 500 in (III.64i are associated with the string correlation length at decoupling 
(perturbations with i > 500 can only be causally seeded after last scattering). Here, the non-separable nature and very 
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FIG. 8 The reduced CMB bispectra for several non-Gaussian models, including (top panels, left to right) equilateral, local, flattened 
models and (bottom panels) warm, feature, cosmic string models (see main text). All five primordial models are normalised relative 
to the constant solution ( |III.48[ > and are taken from Fergusson & S hellard| ( |2009| ). The analytic cosmic string bispectrum l |III.64^ is 
multipUed by (^i^2^3)'*'^ and is taken from |Regan & Shellard, ( |2009l [ 

different scaling of the string CMB bispectrum are clear from a comparison with ( |111.50[ ). Moreover, given the late- 
time origin of this signal from string metric perturbations, the modulating effect of acoustic peaks from the transfer 
functions is absent, as is clear from fig. |8] This is just one example of late-time phenomena such as gravitational 
lensing, secondary anisotropics and contaminants which are accessible to analysis using the more general CMB mode 
expansions ( |111.60| l. 
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E. The estimation of /nl from CMB bispectra 

In light of the previous discussion, it is evident how measurements of the bispectrum from CMB experimental 
data-sets are able to provide information about the primordial three point function of the cosmological curvature 
perturbation field at the end of inflation. This in turn allows us to put significant constraints on inflationary models, or 
on alternative models for the generation of cosmological perturbations. We will now start dealing with the problem of 
bispectrum estimation in the CMB, as a test of primordial non-Gaussianity. 

Let us assume that we have measured the three point function of a given CMB dataset. There are now two general 
ways to exploit this information: 

1. Tests of the Gaussian hypothesis. By comparing the measured three point function to its expected distribution 
obtained from Gaussian simulations we can detect if some configurations present a significant deviation from 
Gaussian expectations. The issue with this approach is that it is sensitive not only to primordial non-Gaussianity, 
but also to any other possible source of NG, including those of non-cosmological origin. Original bispectrum 



tests of this kind on COBE maps (Ferreira et al. 1998 1 revealed significant deviations from Gaussianity in the 



data. This NG signature in the three point function seemed to be localized in harmonic space around multipoles 
€ = 16 and was object of much scrutiny {e.g. [Bromley & TegmaAl|1999t[Magueijo||2000a|b[|Magueijo et al.\ 
1999 2000V It was then finally ascertained that the detected signal was not cosmological in origin, but due to a 



systematic artifact (Bandayefa/j 2000). Moreover, the overall statistical significance of the result disappeared 



in a later analyisis involving the measurement of all the bispectrum modes available in the map (Komatsu 



et al. 20021 (only a subset of all the configurations had been studied before). General tests of Gaussianity are 



very useful to identify unexpected effects in the data, and to monitor systematics. However, as long as we are 
interested in a primordial NG signal, it is better to follow the approach of making an ansatz for the bispectrum 
we expect from the theory under study and obtain a quantitative constraint on a given model. This approach is 
outlined in the point that follows. 

2. ff^i estimation. In this case we choose the primordial model that we want to test, characterizing it through its 
bispectrum shape. We then estimate the corresponding amplitude f^^'^'^^ from the data. If the final estimate is 
consistent wih f^^'^'^^ - 0, we conclude that no siginificant detection of the given shape is produced by the data, 
but we still determine important constraints on the allowed range of /jIJ'l Note that ideally we would like to 
do more than just constrain the overall amplitude, and reconstruct the entire shape from the data by measuring 
single configurations of the bispectrum. However, the expected primordial signal is too small to allow the signal 
from a single bispectrum triangle to emerge over the noise. For this reason we study the cumulative signal from 
aU the configurations that are sensitive to f^{^'^^^. 

Since in this review we are concerned with the study of the primordial bispectrum, we will take the latter approach, 
and deal with the problem of /nl estimation from measurements of the bispectrum in CMB maps. We will first present 
a cubic estimator that optimally extracts the /nl information from the data contained in the bispectrum (section [lII.E. 
We will then address the issue of understanding if this optimal cubic statistic extracts all the possible information 
available on /nl in the data, or if there is enough additional information beyond the three point function to allow 



more precise /nl measurements using non bispectrum-based estimators of /nl (section III.E.2 1. We will then discuss 



concrete numerical implementations of bispectrum estimators (section III.F i and review the experimental constraints 



on /nl obtained from bispectrum analysis of WMAP data (section III.G i. Using a standard Fisher matrix analysis, 
forecasts on the /nl error bars are achievable for future CMB surveys (section III.H i. Following, we will study the 
NG signals in the map that could contaminate the primordial NG measurement and how they are dealt with when 
analyzing the data. Finally we will describe algorithms for the simulation of primordial NG CMB maps that are useful 
for testing and vahdation of estimators before applying them to real data. 

In the following, we assume that the reader is familiar with essential concepts in statistical estimation theory, such as 
the definition of a statistical estimator, the role played by maximum likelihood estimators in statistics, the definitions 
of unbiasedness and optimality, and the definition and main applications of the Fisher information matrix. The reader 
unfamiliar with these concepts can consult the appendix of this review and references therein. 
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1 . Bispectrum estimator of /nl 

In this section we are concerned with the statistical inference of /nl from measurements of the bispectrum of the 
CMB anisotropies. We recall that we defined /nl earlier as the amplitude of the bispectrum of the primordial potential. 
In principle, we can include both temperature and polarization multipoles a^^j^ in the analysis, in order to maximize 
the available data. However, for clarity we will consider only temperature multipoles in the following, and omit the 
superscript T in Of,,,, for simplicity of notation. The extension to polarization is conceptually straightforward, and will 
be discussed in a following paragraph. We will start by considering a simple cubic' statistic written in the form; 

■^^^ ^ "z!? X ^wr'"^i"'i'^^2"'2'^^3m3 • (III.66) 

In the previous equation /.jl represents the statistical estimate of /.jl from the data, are the multipoles of the 
observed CMB temperature fluctuations, are some weight functions, and AT is a normalization factor that has 

to be chosen to make the estimator unbiased i.e. to ensure that: 

</nl> = /nl . (111.67) 
We now want to find the weights that provide the best estimator (i.e. the minimum error bar estimator) 



within the class of cubic statistics written in the form (III.661. It is a well known result (see Appendix A) that the 



best unbiased estimator of a parameter from a given data-set is the maximum-likelihood estimator. In order to answer 
our question we then have to write the bispectrum likelihood as a function of the parameter /nl, and maximize with 
respect to /nl- 

In the assumption that the bispectrum configurations are characterized by a Gaussian distribution^, maximizing the 
likelihood is equivalent to minimizing the following 

(f d/nl = 1 T>ohs 

= Zj ~2 ' ^"^-^^^ 

where B'/^/^^^ is the observed angular averaged bispectrum i.e. , by definition: 

nobs _ ( ^1 ^2 ^3 \ obs obs obs /ttt /-q-v 

and (T^ is the bispectrum variance i.e. the af„, six-point function: 

O"^ = {ai^nnat2m2^hmi<^timi^tsmia(^m^) . (III.70) 

We will now make the assumption that we are working in the weak non-Gaussian limit i.e. /nl is small and the 
distribution of aim can be approximated as Gaussian in the calculation of the variance. The implications of this 
approximation will be discussed in greater detail in the following sections; for the moment it will suffice to point out 
that the weak non-Gaussian approximation is generally a good one since most inflationary models predict /nl to be 
small, and because the level of primordial non-Gaussianity is already constrained to be small by WMAP measurements 



(Komatsu et al. 2009 Smith et al. 2009 1. After restricting indeces so that ti < £2 ^ ^3 and {4 < {s < ((,, the six point 



function above can be calculated using Wick's theorem, yielding (Luo 1994| i: 



) = A C, Q2Ce,6lyl6',l6X6Zt ■ (111.71) 



' The estimator is dubbed cubic due to the fact that it contains the third power of the random variable a[,„ . 

^ This is not strictly true, but it is a good approximation. The same approach applies to most cosmological observables 
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In the last formula A is a permutation factor that takes the value of 1 when all ^'s are different, 2 when two ^'s are 
equal and 6 when all ^'s are equal. We can now substitute (III. 71 1 into ( III.68| l, and differentiate with respect to /nl to 
get an explicit expression for the optimal cubic statistic we were looking for; 



/n 
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N 



Z 



\t,,m,] 



^mim2m, , /nl = 1 



N 



- Z 



(III.72) 



(III.73) 



where bt^e^tj. is the reduced bispectrum and is the Gaunt integral defined by equation (111.39 1; TV is the nor- 

malization factor mentioned at the beginning of the paragraph, that guarantees the unbiasedness or the estimator. 

Note that the noise and window function of the experiment are included in the C/ and bc^c^i^ that appear in the 
formula above, with the following replacements: 



(III.74) 



^ is the window function (not to be confused with the weights W) and A^^ is the noise power spectrum (constant for 
uncorrelated white noise). The noise is assumed to be Gaussian, thus characterized by a vanishing three point function. 



Comparing our result (III. 72 1 to the initial ansatz (III. 66 1, we then see that the optimal weights are: 



(III.75) 



in other words we are weighting the observed bispectrum by its expected signal-to-noise ratio. 

We have now constructed a statistic that optimally extracts the information about /nl from the bispectrum of the 
map. The question now is: is there additional information about /nl in the map that is not contained in the bispectrum? 
This issue will be investigated in the following sections. For the impatient reader we anticipate that the answer is no: 
the bispectrum statistic built here is actually the minimum error bar estimator of /nl from CMB data. 



2. Optimality of the cubic estimator 



In this section we address the issue of whether the cubic statistic (III. 72 1 optimally extracts all the /nl information 
contained in the ag,,,, or if other statistical estimators (e.g. four-point function, or pixel space statistics such as the 
Minkowski functionals, or again wavelet estimators, just to make a few among many possible examples) are able to 
produce smaller error bars and are thus more efficient than the bispectrum. 

In a non-Gaussian primordial CMB map, the acm likelihood depends on the NG parameter /jl- We will indicate 
it with p{a\fNL), where a indicates a vector including all the flf,„'s (we will assume all other cosmological parameters 
are fixed, and concentrate on /nl)- It is a well known result in parameter estimation theory that there is a lower limit 
on the error bars that can be assigned to a given parameter (in our case /nl)- Such lower limit, also known as the 
Rao-Cramer bound, is defined in term of the Fisher matrix F as^: 

A/nl > , ^ . (III-76) 
We remind the reader that the Fisher matrix is defined as: 

/52ln/?(a|/NL)\ /TTT'7'7^ 



^ We again refer the reader unfamiliar with these concepts to the brief summary provided in appendix A 
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If we can show that the bispectmm estimator of the previous section saturates the Rao-Cramer bound for the Of,,, Fisher 
matrix above, then we conclude that it provides the best (i.e. minimum variance) estimate of fyi^fmm the data, rather 
than just the best /nl estimate /rom the bispectmm of the data. In other words, no more information about /nl could 
be extracted from the ac,„ than the information contained in the bispectrum. The aim of this section is to show that this 
is actually the case. 

The issue of the optimality of bispectrum estimators of /nl was addressed in great detail in Babich ( 2005[ l. In 
this section we will basically review the main results of that study, referring the reader to the original paper for their 
complete derivation. 

As we mention in appendix [A] there is a sufficient and necessary condition for an estimator fi to saturate the Rao- 
Cramer bound, expressed by formula ( A.5 i. This condition, applied to our case, reads: 



51n/5(a|/NL) 



= F 



NL 



/nl/n 



^(£(a)-/NL) 



(III.78) 



Our aim is to show that the bispectrum statistic ( |III.72 1 satisfies this condition. We then need to start from a compu- 
tation of the full likelihood p(a|/NL) for a general primordial non-Gaussian model. Following Babich (2005) we will 
start by limiting ourselves to the particular case of the local model. 

We recall from a previous section that local NG is the only case for which an explicit expression for the primordial 
potential is provided. In real space: 



cD(x) = O^x) + f°l- [<l>i(x) - (0i(x))] . 



(III.79) 



Starting from this formula it is possible to obtain a likelihood function for O, dependent on the parameter . This 
is done by means of an expansion in terms of the order parameter /jjj°j^ (0^(x)). The full expression for the Probability 



Density Function (PDF) P(<D|/Nf ) (see 



and schematically written as: 



Babich 



2005 I can be expanded around its Gaussian expectation for /J 



■loc. 
NL 



In Pmf^ ) = In PcmC) + f^l- In PNcmC) + O (f^^ (oi(x)>') , 
where C is the covariance matrix of the Gaussian part of the potential O i.e. 



C = <Oz,(xi)(Dax2)) . 



(III.80) 



(III.81) 



Formula pil.80[ ) is then telling us that the logarithm of the full likelihood can be decomposed into the sum of a 
Gaussian likelihood Pq, plus a NG term that depends linearly on f^^-, and that this decomposition is accurate up to 

terms of order O (f^i^ ^<l)^(x)^ j, i.e. we are assuming that NG is weak, as we did in the previous section. 

After computing ^"(01/^°!^ ), one has to account for 2D projection and radiative transfer in order to obtain the required 

in spherical 



Babich 



likelihood P(a|/^°j|J ). As shown in 
harmonics and performing the functional integration: 



( 2005 I, this can be achieved by expanding the PDF 



III.80 



^(al/NL) = J d''<bd^^> atm- J drr^at(r)<:>(Jr) 



(III.82) 



where 6^^ is the Dirac delta function of dimension M, and M < N due to the 2D projection"* The previous formula 
can be derived by recalling equation ( |III.31| l, together with the well known formula in probability theory: 

P(y) = J dxP(x)6D (y - F(x)) , (III.83) 



* As noted in Babich 2005 the additional degrees of freedom do not affect the CMB anisotropics and can therefore be integrated out. 
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where 5o is again the Dirac delta function, x and y are random variables linked by the functional relation y = F(x), 



and P(x), P(y) are the PDFs of x and y respectively. Solving the functional integral (III.82 1 yields (Babich 2005 1: 



lnP(a|/; 
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+ /2(a, /nl) + O (/^L ■ (111.84) 



In the previous formula we can recognize the standard atm PDF, valid in the standard Gaussian case, in the first term on 
the r.h.s. Added to this we find a first order /NL-correction proportional to the CMB angular bispectrum. Higher order 
correlators are not present at order (9(/nl ('1'/,(x))). For reasons that will become clear shortly, although we have not 
computed it, we have explicitly denoted the 0{j^^^ ^<l>^(x)^) term in the expansion with h{a,f^\^). Note that, besides 
assuming weak NG in this formula, we are also assuming rotational invariance (this is evident from the fact that the 



ai,„ covariance matrix appearing in the Gaussian piece of (111.84 1 is diagonal and equal to Ce). Rotational invariance 
is a general property of the CMB sky but it is broken when we deal with real CMB measurement characterized by 
inhomogeneous noise patterns and sky cuts. We will investigate these effects in the following section. For the moment 
we consider the purely ideal case described by (III.84i. Armed with the PDF expression for local NG, and recalling 



the necessary and sufficient condition (III.78 i, we can finally determine whether the estimator (III. 72 1 is optimal or 
not. First of all we see that: 



d\np{a\f°n 
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■0{fl^{^l{^))). (III.85) 



We then see from combining (III. 85 1 and (III.72 1 that: 



glnj.(a|/^°i;) 
<L 



/NL(fl) + 



^/2(a,/NL) 

5/nl 



(III.86) 



We now see that, in order for the necessary and sufficient condition for optimality ( III.78 1 to be verified, we need the 
(dh/dfyi,) term to be exactly equal to -/nl- The second order quantity I2 should then be calculated explicitly in the 



expansion (III.84i in order to complete the calculation and verify if, or under which conditions, this is true. However 
this turns out not to be necessary if we consider the following "regularity condition for a PDF".^ For a general PDF of 
a random variable x depending on a parameter A we have 



din p(x\A) 
dA 



= 



(III.87) 



Since this regularity condition must be valid for each value of the parameter A (A ^ f^^- in our case), it is clear 
that it must hold tenn-by-tenn, i.e. at each order, in the expansion (III.84i. By taking the average value of equation 



(III. 86 1, keeping in mind that the estimator is unbiased, and imposing (III. 87 1, we then find that the average value of 



c?/2(a, /.jl)/5/nl must be exactly equal to -/nl- If we could then replace i9/2(a,/NL)/<5/NL in (III.86 1 with its average 
value we would exactly obtain the condition for optimality and conclude that the cubic estimator (III. 72 1 saturates the 
Rao-Cramer bound. For present CMB experiments, the terms in the expansion (III.84i are evaluated summing over 
a large number of /"-modes^, or equivalently in pixel space, averaging over a large number of pixels (~ 10^ and 10^ 
for WMAP and Planck respectively). For this reason we expect that the error made by replacing the /nl -order term 



in (III. 86 1 with its average value will be very small. In (Babich 2005 1, an estimate of this error has been done in the 



' The condition jlll.87[ can be easily derived remembering that, for a given random variable x with probability density p(x), we have by def inition 



(f (x)) = J dxF{x)p{x), and substituting F(x) — » 51n p{x\A)/dA in the previous expression. One then finds that the regularity condition III. 87 
holds, provided the order of integration and differentiation can be exchanged (hence the "regularity condition" qualification). 
* (max - 500 for WMAP, {^ax - 2000 for Planck in the signal dominated regime. 
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approximation of neglecting radiative transfer and projection effects (i.e. working in 3D with the primordial potential, 
rather than with a CMB map). The conclusion was that for a number of observations > 30 the approximation above 
works very well. Moreover the variance of the /NL-order term scales like 1 /N. In the full radiative transfer case we 
expect the scaling to be unchanged, although the coefficients in front of it that led to the > 30 estimate might change. 
However, as noted above, the number of pixels in present-day experiments is many orders of magnitude larger than 30. 



That leads us to conclude that the approximation of replacing the average of the first order term in equation (III.86 1 is 
a very good one. We then reach the following important conclusion: 

for a rotational invariant CMB sky, in the limit of weak NG, the cubic estimator defined by formula ( |III.72| l 



is the best unbiased CMB estimator of f^^' 

Let us now move to problem of generalizing the last conclusion to shapes different from local. In this case a full 
expression of the primordial potential <l)(x) is not available. The steps that lead to the conclusion that the local /nl 
estimator is optimal can thus not be reproduced. However it was pointed out in |Babich| ( |2005| l that, in the limit of 
weak NG, the full CMB NG likelihood can still be expressed in terms of its power spectrum and bispectrum by mean 
of an Edgeworth expansion, regardless of its full expression. The Edgeworth expansion is basically a way to express a 



NG PDF as a series expansion around its Gaussian part ( Bernardeau ef aZ. 2002 Bernardeau & Kofman 1995 Taylor 



& Watts 2001| ). For CMB anisotropics one finds, at the end of the calculation: 
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(III.88) 



It is easy to see that InF(a|/NL) takes the same form as in equation (III. 84 1. For this reason all the previous derivation 
applies also to the present case and the following conclusion holds: 

in the weak non-Gaussian limit and assuming rotational invariance of the CMB sky, the cubic estimator 
pil.72| l is the best unbiased CMB estimator of fm^ for any non-Gaussian shape. 

Before concluding this section, we would like to stress that, despite the technical complications arising in the 
detailed probe of the bispectrum estimator's optimality, the physical reason behind this result is quite clear. We can 
always expand the a(„, PDF in series of its momenta. The order parameter of this expansion is (/ml (O)). This 
parameter is the natural measure of the amplitude of primordial NG, and it is actually predicted by inflation to be very 
small. For this reason higher order momenta in the primordial non-Gaussian ac„, PDF are suppressed with respect to 
the bispectrum. Basically, the information on /nl in a CMB map is entirely contained in the three point function. 



3. Breaking rotational invariance 

So far the assumption of rotational invariance of the CMB sky has been made when probing the optimality of the 



cubic estimator (III.72i. In an ideal situation, the CMB sky is clearly rotationally invariant. However, two elements 
break rotational invariance in a CMB map derived from a real experiment: an anisotropic distribution of noise in 
pixel space, and a galactic mask. Anisotropic noise comes from the fact that the CMB sky is generally scanned in 
a non-uniform way: regions that are less contaminated by astrophysical foreground emission are generally observed 
more times, and are thus characterized by a lower noise level (see fig. |9] for an example). A sky cut has also to be 
introduced in order to remove the regions on the galactic plane that are most contaminated by foregrounds. When 
rotational invariance is broken the considerations of the previous two sections do not strictly apply anymore and the 



estimator ( III.72| i become sub-optimal. However, the same Edgeworth expansion approach that was adopted in the 



previous section can still be applied, but this time keeping rotation invariance breaking terms in the calculation, in 
order to find the new more general form of the optimal estimator The general estimator turns out to be the sum of two 
terms: the first term is cubic in acn and is analogous to the one appearing in the rotationally invariant case, while the 
second term is linear in a(m and accounts for breaking of rotational invariance. The explicit expression of this general 
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optimal /nl estimator is (Creminelli et al. 2006| l: 

-3 <af,«„a^2m2«^3m3> C^/,„| ^.,,„^C7^'„,^ ^^,„^af„„,) (III.89) 
^ = ^ {"■Cmiahmiaeym^)Ccl„^^f^^ Cil,^^^_ f^„^C^^^^^^^ (III.90) 

[Cm,] 

In the rotationally invariant case the a(,„ covariance matrix Cctmuiinn is diagonal and equal to Cc, while the linear 
term is proportional to a monopole. We then recover the form of the cubic estimator III .12 as expected. Note that in 
the signal dominated regime of the experiment under study {e.g. £ < 300 for WMAP and d < 1000 for Planck), and 
if the mask is not too large, then the simple cubic estimator ( III.72[i is still basically op timal, since we are in a nearly 
rotationally invariant case. For small masks it has been shown in Komatsu & Spergel| (|2001 ) that the bispectrum and 
power spectrum of the map are, to a good approximation, just rescaled by a factor /jj.,,, representing the fraction of the 
sky left free by the mask, i.e. 



unask 



Cmask _ r r^J 



,fuU sky 



(III.91) 



In this case one can then assume the covariance matrix to be diagonal, and account for the effects of the mask by 
correctly rescaling the normalization term in order to keep the estimator unbiased. This nearly rotationally invariant 
estimator then takes the form: 



NL 



{t„m,] 



/^mim2m3 7 /NL = n 

.V y^^itih "uhhl 
^ - J sky 2_i — r r r 



(III.92) 
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4. Large /nl regime 



The approximation of weak non-Gaussianity is the basis for all the results derived so far. One can then ask at 
which point {i.e. for which values of /nl) this approximation breaks down. As we observed earlier, the Edgeworth 



expansion (III. 88 1 shows that the likelihood of a generic primordial NG distribution can be expanded in series of 
its momenta, with order parameter (9(/nl (•^^(x)^). We know that ~ 10"^, while WMAP observations already 
constrain /nl ^ 100. That means that the order parameter of the PDF expansion is ~ 10"^ and thus the weak NG 
approximation seems a very good one in the entire range of allowed and predicted /nl- However a subtle effect has 
been pointed out by ICreminelh et aT] ( |2007[ ), that change the previous conclusions in certain cases. Let us quickly 
summarize their main results. We already saw that, for the angular averaged bispectrum of a Gaussian temperature 
field: 



(III.94) 



We then included this expression for the variance in the weights of the optimal estimator (III. 72 1, and in the normal- 
ization factor N. It is easy to see that in this approximation the variance of the estimator can be predicted as: 



(III.95) 



29 



However the approximation of taking /nl = in the calculation of the estimator variance is not always a good one if 
A/nl is dominated by squeezed configurations^, or more in general by configurations in which one of the /"s is small. 
It turns out that, in these cases, the /NL-dependent corrections to the Gaussian expectation of the bispectrum variance 
become important when /nl gets large enough. This effect increases the variance of the estimator with respect to the 
expectation for /nl - 0. There is a clear physical interpretation for this. One can see, for example by calculating 



(III.95I in the simple Sachs Wolfe case (see also section III.Hi that the variance of the estimator scales roughly like 



l/^mn.v, or equivalently like 1 /A^y„,v in pixel space, Npi^ being the number of pixel in the observed map. This increase of 
the signal-to-noise ratio with the number of pixels is due to the fact that more and more bispectrum configurations are 
included into the sum over modes to estimate /nl- However, if the signal is completely dominated by low-^ modes, as 
in the local case, there is an intrinsic large cosmic variance, due to the small number of low-^ configurations. Clearly, 
cosmic variance cannot be beaten by increasing the resolution of the map. Creminelli et al. ( 2007) then found that 



for Npix getting large, the S IN of the estimator for local NG grows asymptotically as (inNpix) i.e. much slower than 
the expected Npi^, that one would obtain by neglecting /NL-dependent corrections in the calculation of the estimator 
variance. They carried out a calculation of the estimator variance in the flat-sky approximation, and neglected the 
transfer functions, to find the following expression: 



((A£)2) 



1 



1 + 



AN ■ ^ 



nln Nn 



(111.96) 



where A is the bispectrum amplitude. We clearly see from this formula what we were stating above i.e. when /nl 
gets large, the variance starts scaling like 1/lnA^^^.^. The same formula also shows the technical point behind this 

behaviour: in the correction term, the order parameter is actually not /nl^'^^ anymore but rather fm.A^^^Npix. This 
enhancement by a factor Npix can make the first order corrections non-negligible anymore. The natural question is 



now how large an /nl we need to make the correction term important in (III.96 1. Following CreminelU et al. (2007 1 
let's call (To the standard deviation of the estimator computed for f^^- - 0. Let's say that for an experiment at a given 
angular resolution (defined by l,„ax in harmonic space or by Npix in pixel space) a value of f^^- is measured, equal 

formula (??)behavior, and 

following relative correction to ctq: 



to ncTQ. Substituting this value into formula (??)behavior, and calling cPj^^ the real estimator variance, one finds the 



cr 



- 1 = 



2«2 



(III.97) 



For an experiment like WMAP the r.h.s. term becomes ~ 1 when /nl is about ba-Q away from the origin. For 
an experiment like Planck /nl has to be about S.Sctq. The definition of an high-/NL regime is thus dependent on 
the experiment under study, as a consequence of the fact that the enhancement of first order terms in the variance 
expression ((Afi)) depends on Npix. In other words we can conclude the following: 

'//nl' '^''^ detected at several cr (in terms of the Gaussian expectation for the standard deviation), 
then f^i^-dependent correction terms in the estimator variance will have to be taken into account, and the 
simple expansion \nL84\ of the CMB likelihood does not constitute a valid approximation anymore. 



One caveat in all this discussion is that formula (III. 96 1 was obtained in flat sky and neglecting the transfer functions 
and it should be checked how dependent the final results are on these approximations. Since the scaling argument is 
based on a very general physical reason, i.e. the weight of squeezed configurations in the local /nl estimator discussed 
earlier in this section, one expects that the general scaling with Npix obtained in ( III.96| l does not depend on the details 



of radiative transfer and 2D projection. Liguori et al. ( 2007^ actually checked the results of this section numerically. 



by applying an implementation of the optimal cubic estimator pil.72| i to full-sky simulations of CMB local NG maps 



^ We recall that by squeezed configurations we mean triangles in which one of the sides is much smaler than the other two, i.e. Ci « {2, £3. 
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with different A^, 



including the full radiative transfer . Although, as expected, the coefficients in formula 



(III.96I change with respect to the simple flat sky no radiative transfer approximation, the scaling of the error bars 
with Npu- follows very well the expectations, going from ~ IjNpix InNpix at low /nl to 1/ In A^p^ when /ml is detected 
at high significance. Since in the large /nl regime the variance starts to scale very slowly, like 1 / In^ Npix, one is led 
to wonder whether the estimator discussed in the previous sections becomes suboptimal at this point, and whether a 
better one can be found. The answer to this question is not immediate. In order to check for the optimality of an 
estimator, as we have seen, one has to see if it saturates the Rao-Cramer bound. However, also the local /nl likelihood 
and Rao-Cramer bound estimated in the previous sections have to be recomputed, since they were obtained neglecting 
higher order terms in /nl'4^'^^- In order to account for the /nl^'^^A'pii enhanced terms it is necessary to produce an 
exact expression of the full likelihood. This can be extremely challenging in the full radiative transfer case, but it 
is feasible in the flat sky no radiative transfer approach that we are considering (and that we showed earlier to be 
a good approximation as long as scaling arguments are involved). Creminelli et al. (20071 proceed to calculate the 



full likelihood in this approximation and conclude that the optimal cubic estimator of weak local NG indeed does 
not saturate the Rao-Cramer bound in the high-/NL regime. The estimator (III. 72 1 is thus no longer optimal in this 
case. They then proceed (always in the flat sky no transfer function case) to derive a cubic estimator that saturates the 
Rao-Cramer bound also for large /'°j^ . We won't enter into the details of this derivation here, referring the reader to 



Creminelli et al. ( 2007 1 for a complete discussion. The main aim of this section was to show under which conditions 
the optimality of the cubic estimator that we discussed in pr evious sections is valid. Since current bispectrum analysis 



of WMAP data (Komatsu et al. 



2009 



and the cubic estimator we derived ear 



Smith et al. 2009) find /^™ ~ 2cro, the weak NG approximation applies, 
ier is indeed an optimal estimator in this case. However, if future Planck 
measurement will produce a detection of /^°j^' at high significance the estimator will have to be modified in order to 
account for the enhanced variance in the high-/NL regime. This is not necessarily a remote possibility if one considers 
that the present central value of /nl from WMAP estimates would produce a ~ Sctq detection with Planck. Before 
concluding this section we would like to remark once again that the variance enhancement discussed here applies only 
to non-Gaussianity of the local type, whose bispectrum is dominated by squeezed configurations, affected by a large 



cosmic variance. For the other shapes the cubic estimator ( III. 72 1 is optimal in both the small and high-/NL regimes 



F. Numerical implementation of the bispectrum estimator 

In this section we turn to the problem of finding an efficient numerical implementation of the optimal bispectrum 



estimator (III. 89 1. Let us repeat its expression here for convenience: 

\t„mi] 

-3 {af,mifl<'2m2«^3m3> Q/,„i/,„2Cfj„3,4m,af4m4] ' (III.98) 
^ - {atimxahmi'^him) QiLi.f4»i4^f2m2/5m5^f3m3,f6m6 {O-Umfihrnsahm,^ ■ (III.99) 



If ;,»!;) 



We remind the reader that this is the full expression, valid for the general non-rotationally invariant case. For a 
rotationally invariant CMB sky the linear term in the formula above vanishes, and the covariance matrix is diagonal 



For details about the numerical implementation of the optimal cubic estimator, and about the generation of NG CMB maps, see sections |III.F| and 

Imjl 



31 



and reduces to Ce, giving the simplified expression ( |III.72| l that we reproduce again here for convenience: 

-,mjm2m2 . /nl=1 



^^'''> = 77 / , ^ ^' " a(,„,^a(^m^ae,m, , (Ill.lOO) 



2 



/^mim2m3 7 /nl = 1 Y 

In a schematic way, the full estimator can be written as: 



TV = y ^ '''''' ''''''' . (m.ioi) 



fi(fl) = , (ni.l02) 

where the "cubic" term is the one containing the product flf,,„,af,„,,flf3,„3, while the linear term is the one dependent 
on a single ag,,, and vanishing in the rotationally invariant case, where it is proportional to a monopole. It was shown 
before in a formal way that a pure cubic estimator becomes suboptimal when rotational invariance is broken, and 
adding the linear term is necessary to restore optimality. It is useful to try to understand the reason of this effect 
qualitatively and in a more intuitive way. Let's assume that we have a map characterized by non-stationary noise, and 
we are observing a region of the sky that was sampled many times so that the noise level in this area is low. That 
implies that the level of small scale power in this large region is lower than average. Now, for a specific realization of 
the CMB sky, this modulation of small scale power on a large region can look like a non-Gaussian signal sourcing a 
squeezed configuration of the bispectrum. On average this effect must cancel if the underlying noise model is Gaussian. 
However, this "confusion" between signal and noise increases the variance of any estimator of a primordial NG signal 
that is peaked on squeezed configurations. We know that this happens for the local model. This heuristic argument 
thus shows that, even though in principle a linear term must always be included when rotational invariance is broken, 
for a realistic noise model only local non-Gaussian estimates will be affected. 



1 . Primary cubic term for 



Let us focus for the moment on the rotationally invariant case, where the linear term vanish, and the covariance 
matrix is simply C - Cc6c^e^5mtm2- We immediately see that a brute force implementation of equation (III.lOOi, 



consisting in computing and summing over all the bispectrum configurations, would take O (imax) operations, where 
{max, the maximum multipole in the calculation, depends on the resolution of the experiment. As mentioned earlier, 
{max ~ 500 for WMAP and {^ax ~ 2000 for Planck in the signal dominated regime. At these resolutions a brute force 
approach would be absolutely unfeasible for a general shape. If however we assume that the primordial shape under 
study is separable then the dimensionality of the problem can be reduced and the overall number of operations scaled 
down significantly, making the computational cost affordable. Let us illustrate this point more in detail. Substituting 



(111.46 1 into the estimator expression (111.72 1, and remembering the identity (111.39 1, 



^mim2fn2 



,i,t2h ^ danYt,,n,{r^)Yt2n,,{XL)YhmA^) , (ni.l03) 

it is possible to rewrite ( [111.1 00[ ) as follows (or more in general as a linear combination of terms of the following kind): 

,,„,Xf,(r)yf,,„,(n) at^,n2^t^{r)Y i,jn^{VL) ^ ac^^;Ki^(r)Y 



&(d) = ^ r dr? fd^i>Yj — 



2 



Ce. 



Z 



Ce 



■ -I- perm. . (HI. 104) 



From an inspection of previous formula we see how, as a direct consequence of separability, the initial sum over the 
indeces ^i^2^3,'Mi'M2'«3 has been factorized in the product of three sums, each running over two indices i,m. This 
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greatly reduces the computational cost from O {€^„ax) 1° ^ {^max) operations. If we define the new quantities 

Mx(r,n) = 2^ Ffm(n), 



MY(r,h) = ^ 

fm 



Ci 



Mz(r,n) = ) — Y(„,(n) , 

1^ 



Ce 
aemZeir) , 



we can recast the estimator expression above in the following form: 



S(a) 



dQ.nMx(r, n)My(r, n)Mz(f, n) + perms. , 



(III. 105) 



(HI. 106) 



where it is evident that we are now calculating our statistic in position space rather than in pixel space. Note how 
the filtered maps Mx, My, Mz can be efficiently calculated using a fast harmonic transform algorithms such as those 
included in the HEALPix package. This fast position space algorithm was initially introduced in |Komatsu et al.\ 
( 2005[l in the cont ext of local /nl estimation, and applied to the estimation of WMAP 1-year data by the WMAP team 



Komatsu et al. (2003 1. It was then applied to equilateral /nl estimation for the first time inl Creminelli et a/.] ( |2006) l. 



An alternative numerical implementation with respect to the one used by the aforementioned authors was introduced 
in [Smith & Zaldarriaga| ( [2006 1. Although different under many technical aspects, this second algorithm is still based 
on the calculation of the position space statistic (III. 106 1; we refer the reader to the original work for additional details. 



This second implementation has been used to produce alternative estimates of 



and /^^ WMAP data, and to 



estimate the amplitude of the orthogonal shape, recently introduced in Senatore et al. ( 2009) . 

Let us now discuss the possible limitations of this numerical approach. As noted in section [llj the separability 
condition is in principle quite restrictive: the only separable shape arising directly from primordial models of inflation 
is the local one. On the other hand, it is still possible to study non-separable models by finding separable shapes that 



are highly correlated to the primordial one. As observed in CremineUi et al. (2006i; Fergusson & SheUard (2009); 
Smith & Zaldarriaga ( 2006| l, the /nl limits obtained from a highly correlated separable shape in this way will be very 



close to those that would have been obtained using the original non-separable model (see again sections II III.B and 



II. B for a detailed discussion of this issue). We know from earlier sections that the other two shapes mentioned so 



far in this section besides local, namely the equilateral and orthogonal shape, have actually been derived as separable 
approximations of theoretical inflationary shapes. These approximations were obtained in an heuristic way i.e. an 
educated guess of a good separable approximation of the shape under study was made, and the correlation was checked 
a posteriori. There is obviously no a priori guarantee that this approach would be easily repeatable for all the shapes 
of interest. The eigenmode expansion method introduced in ( [Fergusson et al.\ [2009| l, and summarized by equation 
( II.26| l, however, provides a general and rigorous method to find separable approximations of any shape, thus enabling 
the estimation of any possible primordial model. In this case, recall that we expand our (non-separable) primordial 



shape function in terms of the separable basis functions Q„ (see ( 11.26 1, constructed from symmetric polynomical 
products qp{k), as 

S{kx,k2,h)^^aprsqp(ki)qr(k2)qsih), — > bi^u^ ^ A\f^i^^ap,., ^ ^dxq^^'^qf^ql^ , (III.107) 



where the second expression for the reduced bispectrum /jf,f,;, ( III.56| l expands in convolved basis functions (III. 57 1 
in harmonic space with 

q^ix) = ^ J dkqp(k) Ai(k) j,(kx) . (III.108) 
In the mode expansion approach, then, the /NL-estimator for a specific model generalises to the following 

&(a) = ^ "P" J '^'"'"^ / ^^P^''' ^'l'^'"' ' (HI. 109) 
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where the fihered maps or shells Mp(r, n) are defined by 



Defining the integral fiprs = J drr^ J dQ.fi M^p Mr Mj), the estimator collapses into the compact form 



(HI. 110) 



prs 



(III. Ill) 



and 



where it is possible to show a precise relationship between the theoretical bispectmm expansion coefficients Up^s 
expectations for the observed coefficients fiprs- 

It was also pointed out in Fergusson et al. (2009) (see also Fergusson & Shellard (2007 1), and summarized in 



formula ( III.60| l, that the separation can be performed directly in harmonic space on the reduced bispectmm bg^^jj, 
rather than on the primordial shape S{ki,k2,k^). This provides an alternative, but equivalent, late-time /nl -estimation 
pipeline with respect to the primordial shape separation approach given above ( III. Ill i. In fact, since orthonormality is 



more direct on the harmonic domain without the intervention of transfer functions, the approach is considerably more 
straightforward conceptually. In this ca se, expectations for the observational expansion coefficients in the orthonormal 
frame 'R„ (with n «-> [prs] see (III.6O1) become simply (ySJJ) = o-JJ, that is, for an ensemble of maps possessing the 
theoretical bispectrum described by the coefficients aj^. This means that for a NG bispectmm signal of sufficient 
significan ce we can consider directly and efficiently reconstructing the bispectrum from the observed coefficients 
ySJJ using (III. 60 1. We also note that the harmonic space separation scheme also allows for the estimation of non- 



inflationary late-time bispectra, such as the bispectrum of cosmic strings, as well as other secondary anisotropics. 



We can then conclude, in light of these developments, that the fast cubic statistic ( III. 1 06 1 can be applied in complete 
generality to any model of primordial NG, as well as to any other potential source of CMB NG. We also point out 
that alternative approaches have been considered for harmonic space /nl analysis using wavelets and binning. For 
example, Bucher et al. (|2009jl recently proposed using a suitable binning scheme in which the full expression for f^f, 



is calculated in a subset of all the triples £i,{2, £3, small enough to make the calculation feasible while maintaining 
calculation accuracy. Approaches based on a harmonic space separation scheme, of course, require the full calculation 
of the reduced bispectrum bg^f^e^ in order to determine the correlation between the theoretical prediction and the final 
expanded or binned bispectrum. The calculation of bf^i^j^ implies the necessity of numerically solving the radiative 
transfer integral (III.30I for all the configurations ^i,/'2,^3 which appears to be intractable in the non separable case, 
since the dimensionality of the problem cannot be reduced. However, this can be achieved efficiently in the general 



case either using the separable mode expansion integral (III.56i or else the hierarchical adaptive approach ^Fergusson] 
|& Shellard] ( |2007l ) discussed in section jllLDl 



2. Linear correction term for /nl 

Let us now consider the realistic situation in which inhomogeneous noise and a sky-cut break rotational invariance 
(see fig.|9]l. In this case two complications arise: 

1. A linear term in Of,,, has to be added 

2. The flf,„ covariance matrix is now no longer diagonal. The inverse covariance weighting C~^a that appears in 



expression (III. 98 1 is hard to compute numerically for high angular resolution experiment, since its size makes 



a brate force numerical inversion impossible. 



A first approach, introduced in Creminelli et al. ( 2006 1, is to simplify the problem by assuming that the covariance 
matrix is diagonal in the cubic term of the estimator, and then finding the linear term that minimize the variance under 
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FIG. 9 Left panel: the KQ75 galactic and point source mask used for non-Gaussian analysis of WMAP data. Right panel: 
anisotropic distribution of the noise for a coadded map of the WMAP V and W frequency channels. These two features break 
rotational invariance of the observed CMB sky, and spoil the optimality of the standard cubic statistic unless an additional 

"linear term" is included, as explained in the main text. Data from the LAMBDA website, http://lambda.gsfc.nasa.gov/index.cfm_ 



this assumption. In other words, we keep the cubic term in the form (Ill.lOOl and compute the variance of this term, 
relaxing the assumption of isotropy at this point'^ 

It turns out that the variance is minimized (while leaving the estimator unbiased) for the following choice of the 
linear term: 



Un 



3_ 

'n 



^mim2m3 I /nl=1 



(III. 11 2) 



where N = fskv Efi^f, C c^C c^C is the normalization term. Despite being sub-optimal with respect to a full 

implementation of equation ( III. 98 1, this choice of linear term has been shown to significantly improve the error bars 



with respect to the simple cubic statistic ( |III. lOO i in presence of a nisotro pic noise. At the same time, the simplicity 
of this implementation in comparison to the full optimal statistic (III. 98 i is manifest, since no C"' terms appear in 



equation (III.112i. Let us consider again a separable primordial bispectrum shape that can be written S {k\,k2,kT,) - 
X{ki)Y(k2)Z(ki) + perm.. Applying the same procedure as we did for the cubic term, the linear term can be recast in 
the form: 



^ i,^t\m\ ^iim^) ^^3/713 

-I- perm.) , 

where we explicitly wrote Cf,,„,/,„,j as {a(^,„^a(^j„^) . This last formula can be rewritten as: 



(m.ll3) 



JUn 



'n 



^ ^ [m 



Zfir) 

+ ^{Mx{r,n)MY{r,n)) 



^'^'^ {Myir, n)Mzir, n)) + ^ {Mxir, n)Mzir, n)) + 



Cf 



(m.ii4) 



' This means that, when we apply the Wick's theorem to the Of,,, six point function in the calculation of the cubic term variance, we take 

{aeim^ae^„^^ = Cl^m^/^mj instead of ^afjmjafjm^) = Qimi ^«if2<5niim2 
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FIG. 10 Error bars (obtained from Gaussian simulations) for the pure cubic (triangles) and pseudo-optimal (stars) implementations 
of the bispectrum estimator, to be compared to the solid red line, representing the Fisher matrix (Rao-Cramer) bound, saturated 
by the full optimal statistic described in the text. Left panel: error bars as a function of the maximum multipole included in the 
analysis. Right panel: error bars as a function of the fraction of the sky considered in the analysis. This analysis included both 
temperature and polarization data. From |Yadav et al.\\200S\ . 



Like for the cubic part of the estimator, we have rewritten the linear term as a fast position space integral. The 
ensemble averages appearing in the last formula can be computed as Monte Carlo averages over a large number of 
Gaussian realization of the CMB sky, characterized by the same beam, mask and noise properties as the experiment 
under study. This pseudo-optimal, but relatively straightforward, imp lementation of the linear term has been adopted 



by a number of groups in order to estimate f^,'^' from WMAP data (Creminelli et at. 



2006 



2007 



Komatsu et al. 



2009 ; Yadav <& Wandelt , '2008'). The full optimal estimator (III.89 1 was implemented only quite recently in Smith] 



et g/. ( 2009 1, where the auth ors developed an effi cient conjugate gradient inversion (e.g. Shewchuk 1994 ) algorithm 
based on earlier results from Smith et al. (2007 1, in order to compute the C 'a pre-filtering in reasonable CPU-time. 
Note that after the inverse covariance matrix pre-filtering is calculated, the numerical implementation of the estimator 
is very similar to the one outlined above for the pseudo-optimal case. The new position space statistic is obtained from 



formulae (III. 106 1, (III.l 14 1, by making the following replacements, wherever the corresponding quantitites appear: 



a 



filtered 

{m ^ 



Mx(r, ft) ^ Mx(r, n) = ^ ac„Xe(r)Yi„,(n) 



(HI. 115) 



with analogous substitutions to be made for the Y, J/, Z, X terms appearing in the same equations. The improvement in 



error bars from the pure cubic sub-optimal estimator, to the pseudo-optimal and optimal statistics is shown on fig. 10 



G. Experimental constraints on /nl 



In order to obtain an estimate of /nl from a given data-set one has first to generate sets of Gaussian CMB maps 
and obtain the MC averages that appear in the linear term expression ( III. 1 14 1, after an inverse covariance pre-filtering 
of the full optimal estimator is implemented. The normalization term N can be pre-computed using formula ( III.46[ ) 
to evaluate numerically the theoretical bispectrum shape for the model we want to estimate. The statistic ( III.98| l can 
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then be computed for the experimental data Oobs to get our result: 



N 



The error bars are then obtained by running the estimator on simulated Gaussian maps 



10. 



jm 



[fNLiOsim)) 



(HI. 11 6) 



(III. 117) 



where (.)mc indicates the MC average and aj,„, a vector of simulated multipoles (obviously including mask, beam 
and noise features of the experiment). For an accurate step-by-step description of an /NL-analysis of WMAP data, 
including details about channel coadding, noise model, beams and pixel weighting schemes, we refer the reader to 
the explanations contained in Komatsu et al. ( 2009| l. The most stringent limits so far have been obtained by applying 
the bispectrum estimator to the WMAP datasets. Constraints have been put on the local, equilateral and orthogonal 
shapes. The best constraints come from the full implementation of the optimal estimator done in Smith et al. (2009 1 
and|Senatore et al.\^m9\. They ai-e, at 95% C.L., 
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(III. 11 8) 
(HI. 119) 
(HI. 120) 

Since the first release of WMAP data, different groups have used the cubic statistic described in the previous paragraph. 



either in its pure cubic form (III. 106 1, or in the improved version including the pseudo-optimal linear term implemen- 



tation (III.l 14 1. The results of different analysis of the WMAP 1-year, 3-year and 5-year datasets are summarized and 



commented in table III where just the local and equilateral shapes have been included since the only constraint on the 



orthogonal shape to date has been already mentioned in ( III. 1 18 1 



H. Fisher matrix forecasts 

The fisher matrix, defined as the curvature of the likelihood function calculated in its peak reassessment (see equa- 
tion (A. 3 I in Appendix [A]), plays a very important and well-known role in parameter estimation theory, not only 
because it defines the optimality of estimators through the Rao-Cramer bound, but also because it allows us to esti- 
mate a priori what the smallest error bars attainable will be for a given parameter (see again Appendix [A). In other 
words, using the Fisher matrix we can forecast how well a parameter will be measured by a given experiment. This is 
very useful in order to optimize the experimental design to the detection of the parameters of interest. In our specific 
case, a Fisher matrix analysis will help us to understand what is the smallest /nl detectable in principle using different 
CMB datasets, and which experimental features can be improved in order to increase the sensitivity to /nl- 



1 . A general derivation 
Formula ( |A.12| l from appendix|A] when applied to our case yields: 

/'R/NL=n2 



/nl/ni 



z 

^1 ^2 ^3 =2 



(III. 121) 



The error bars can be obtained from Gaussian simulations as long as the weak NG approximation applies. As we saw earlier, this works at any 

/nl -dependent. In this case the error bars would need to be 
has been reported, so working with G maps is at this stage 



rloc 
'NL 

calculated from NG simulations of . So far no high-significance detection of f^^^ 



/nl for any shape, except for the local shape when a large f^^- makes the error bars /nl -dependent. In this case the error bars would need to be 



loc. 



sufficient to get accurate eiTor bars. 
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Local 



Equilateral 



rule cuuic 


-58 < /nl < 134 (Komatsu ef a/.|l2003' Wl) 
-54 < /nl < 114 ( Spergel et al. . 2007 , W3) 


-366 < /nl < 238 (Creminelli e/ a/.ll2006l, Wl) 
-256 < /nl < 332 (Cremmelli et al. 2006 , W3) 


Pseudo-optimal 


-27 < /nl < 121 (CreminelliefaZ. , 2006 Wl) 
-36 < /nl < 100 ( Creminelli et al. , 2006 W3) 
27 < /nl < 147 (jYadav & Wandelt 2008 , W3) 
9 < /nl < 129 fSmith et al. . 2009 , W3) 
-9 < /nl < 111 (Komatsu al. . 2009, W5) 


-151 < /nl < 253 (jKomatsu ef a/.||2009( W5) 


Optimal 


12 < /nl < 104 (Smith et al. 


2009 W3) 


-125 < /nl < 435 iSmitheffl/. 2009 W5) 


-4 < /nl < 80 (Smith et aC 


2009 W5) 





TABLE 111 Constraints on f{^["' ,f^'l^''' , obtained by dilFerent groups on the one-year (Wl), three-year (W3), and five-year (W5) 
WMAP data releases. Different rows correspond to the different implementations of the /nl estimator described in the text: the 
"pure cubic" implementation jlll. 100^ in which no linear term is included, the "pseudo-optimal" implementation ( |II1.1 14| l in which 
a linear term is added but the covariance matrix is assumed diagonal in the cubic term, and the fully "optimal" implementation 
|nr98j. As we noted in the text the linear term is important mostly for estimates of local NG, since anisotropic noise "mimic" 
squeezed configuration. For this reason "pure cubic" estimates of equilateral NG in the table above are nearly optimal, while local 
ones are significantly sub-optimal, especially because they have to be confined to the pure signal dominated region I < 300, where 
the assumption of rotational invariance is correct. There is a certai n degree of friction betwe en some of the results shown. In 



particular the 27 < f^^- < 147 WMAP 3-year estimate obtained in Yadav & Wandelt (2008 i, corresponding to a " nearly 3-cr' 
detection of local NG, seems not to agree well with the 9 < /JjJ^ < 129, ~ 2.3cr result obtained on the same data set in Smith et al.\ 
p009 ). The origin of the discrepancy is unclear, although it is argued in [Smith et fl/.| lj2009) that it might be due to dilterences in 
the coadding scheme of different data channels, or analogous differences in the choice of some weights. As pointed out in SmiTh] 
|ef fl/.|p009) , one additional advantage of the fully optimal implementation of the estimator is actually that all the ambiguity related 
to the use of different coadding schemes disappears, since the optimal coadding strategy is automatically selected in the inverse 
covariance filtering process. Another discrepancy is that between the two equilateral constraints on WMAP 5-year data. It seems 
that the pseudo-optimal estimator produces better constraints than the optimal one. This is clearly not possible. [Smith et fl/.| ( |2009} 
claim that their numerical pipeline calculates the theoretical ansatz for the bispectrum shape more accurately than it was done 
before. That is due to a subtlety that went unnoticed in previous works, consisting in the necessity to extend above the horizon the 
upper integration limit in the calculation of the equilateral shape related quantities y3f(r), ydr), S[{r) (see equation (??)). This is 
required in order to obtain stable numerical solutions, and it calls for a reasessment of the expected and measured error bars, that 
actually increase with respect to previous calculations. 



where 5^,6^, is the angular averaged bispectrum {i.e. the measured quantity). This can be rewritten in terms of the 
reduced bispectrum as: 



/nl/nl 



where we have defined (see also w/,/,/, in (III. 58 i): 



1 y /2 liiMiL, (ni.i22) 



(2^1 + 1)(2^2 + 1)(2^3 + 1) Ui d (. 



h,.j, - y " - • ' \ ) ■ ("^-123) 

Note how the features of the experiment enter the Fisher matrix through the parameter tmaxi defining the angular 
resolution, and in the angular power spectrum expression in the denominator, that contains the angular beam and 
experimental noise: 

Ce^C[Wc + Nt , (in. 124) 

where Ce is the theoretical power spectrum for a given set of cosmological parameters, Wc is the beam of the experi- 
ment, and Nc is the experimental noise. A^^ is a constant for uncorrelated noise. Likewise, the theoretical bispectrum 
will be convolved by the experimental beam. 

Bi,t,t, = St,e,t,Wt,Wi,Wt, . (ni.l25) 
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Note that, since the noise is generally Gaussian, its three point function vanishes. The experimental noise thus only 
enters in the denominator of the Fisher matrix expression. The effects of partial sky coverage can be easily accounted 
for. From (III.91 1 it follows that if only a fraction f^ky of the full sky is covered then the Fisher matrix takes a f^ky 
factor in front, that produces a degradation of the error bars of ^jfskx- 

We saw previously that for separable shapes the reduced bispectrum can be calculated either analytically, under 
some sir nplifyin g assumptions on the transfer functions {e.g. the Sachs-Wolfe approximation), or numerically through 
formula ( III.43 i. It is then possible to evaluate numerically the fisher matrix and the corres ponding error A /nl = VI IF . 



In the context of /nl estimation, the first calculation of this kind was done for /'°'^- in Komatsu & Spergel (2001 



where it was found that WMAP could reach a sensitivity A/nl - 20 (note how this bound is actually saturated by 
the optimal estimator results presented in table III i, while Planck (The Planck Collaboration 2006 1 could go down 
to A/nl - 5". What allows Planck to improve on WMAP is that it has a much better angular resolution and that 
it is cosmic variance dominated in a very large range of scales i.e. the power spectrum signal C(bt is larger than the 
noise N( up to 4ma - 2000. Angular resolution and sensitivity are the two factors that increase the ability of a CMB 
experiment to constrain /.jl. This information is provided by the Fisher matrix expression (III. 122 1. Looking at such 



expression, we notice how the signal-to-noise ratio is obtained by adding over all the bispectrum configurations up 
to (max, weighted by their variance. Thus, the higher ^,„f„ is, the more configurations are included in the sum and 
the larger is the final sensitivity to /nl- On the other hand we see that, if the power spectrum of the instrumental 
noise appearing in the variance term in the denominator dominates from a certain Cs=n, then the signal contribution 
is suppressed above that threshold by the noise power spectra appearing in the denominator of (III. 122 1. So what 
determines the sensitivity of a CMB experiment to /nl is the range of ( over which the instrumental noise is low, so 
the experiment is cosmic variance dominated. This range is ^ < 2000 for Planck and { < 500 for WMAP, hence 
Planck can obtain tighter constraints than WMAP. This is shown in fig.[TT[ where Fisher matrix forecasts of /.jl are 
plotted for different CMB experiments: the predicted error bars decrease with t up to the angular scale at which the 
measurements start to be noise dominated, after which the /ml signal-to-noise ratio saturates. A simple calculation 



done in Babich & Zaidarriaga ( 2004,) taking the Sachs Wolfe approximation, and working in flat sky, showed that 
before noise dominates the signal-to-noise ratio for the local shape grows as: 



In 



(HI. 126) 



where the (In) is dictated by the coupling between large and small scales introduced by squeezed configurations, from 
which most of the local signal comes. 

Note also how, in absence of experimental noise, the beams in the numerator and in the denominator of ( |ni.l22[ ) 
cancel each other out. An ideal noiseless CMB experiment would then have a signal to noise ratio indefinitely growing. 
However, this would not imply infinite sensitivity to /nl, because above a certain {^ax secondary anisotropies would 
start to dominate. Fisher matrix analysis of the equilateral shape (Sefusatti et al. 2009 Smith & Zaidarriaga 2006 



e.g. ) showed that the minimum achievable error bars in this case are A/nl - 100 and A/,jl 
Planck respectively'^. Additional shapes are studied in Smith & Zaidarriaga (2006jl. 



60, for WMAP and 



2. Polarization 



Babich & Zaidarriaga ( 2004| l showed with a Fisher matrix analysis that the CMB E-mode polarization measurements 



can be used to improve the sensitivity to /nl- Although we have dealt so far only with temperature bispectra and related 
estimators, including polarization is fairly straightforward. As usual the calculation starts from the formula (III. 30 1 



Note that all the errors quoted in this section are at 1 -o". 

Note how the larger error bars in this case with respect to the local constraints does not reflect a higher sensitivity of CMB measurement to f^', 
but only the conventional choice of the normalization of the bispectrum amplitude in the definition of /nl- The normalizations are in fact chosen 
in such a way that the bispectra have the same value for equilateral configurations l\ = (2 = (3, where the local bispectrum is suppressed, and 
the equilateral bispectrum is peaked. 
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linking the multipoles of CMB anisotropics to the primordial potential (5>, but this time including the polarization 
radiative transfer A?(fe) in the convolution integral: 



-0' f ^Af 



W<l)(k)Fak)., 



(HI. 127) 



The bispectrum is then defined in the usual way, but this time more configurations can be built by correlating temper- 
ature and polarization multipoles; 
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(HI. 128) 



The point to emphasize is that the polarization signal is generated on scales where the temperature signal is suppressed 
by Silk damping. The polarization bispectra thus open a window over a new k-range in the 3D — » 2D projection 
k — > and increase the overall information available. In other words, since the new configurations TTE, TEE, etc. 
including polarization are partially independent of the pure temperature (TTT) bispectrum, adding those additional 
configurations to the Fisher matrix (and to the actual /nl estimation from data) increases the total signal available. The 
Fisher matrix expression now becomes: 



pqr ijk tihh 



pqr\ijk Ci ^2^-. 



(ni.i29) 



where i, j, k, p, q and r run over the T and E superscripts. We still work in the assumption that all the quantities involved 
are Gaussian, but now the different bispectra of temperature and polarization are correlated for a given configuration 
^\,{2, £i, thus defining a multivariate Gaussian distribution. The full covariance matrix between bispectra (indicated 
by cov in the formula above) has then to be evaluated. A numerical evaluation of 129)1 shows (|Babich & Zaldar- 



|riaga| [2007 ) that, for an ideal (i.e. noiseless) experiment, adding the polarization signal produces an improvement of 
a factor ~ 2 on /nl constraints. For WMAP, adding polarization bispectra produces very little improvement, since 
polarization data are mostly noise dominated. For Planck, however, including polarization does generate a significant 
improvement, bringing the forecasted error bars from A/nl - 5 to A/nl - 3.5 . Some error bar forecasts from temper- 
ature and polarization bispectra as a function of i,„ax for different experimental designs including WMAP and Planck 
are shown in fig. 1 1 Motivated by this analysis, Yadav efoT] ([2007 2008) have implemented a bispectrum estimator 
of /nl including both temperature and polarization bispectra. All the general considerations about optimality and the 
numerical implementation techniques described in previous sections apply in an analogous way to the temperature + 
polarization case, although the presence of additional bispectra with a non-trivial covariance matrix introduces a few 



additional technical complications. We refer the reader to , Yadav et a/.| ( |2007{|2008j ) for further discussion. 



I. Non-Gaussian contaminants 



So far we have considered only primordial non-Gaussianity as a source of the three-point function of the CMB. 
However many other astrophysical and cosmological effects can produce an observable angular bispectrum. Among 
these, diffuse astrophysical foreground emission (see e.g. Bennett et al. 2003 [ Leach et al. 



2008 and references 



therein), unresolved point sources (e.g. |Komatsu et al.f |2003| l and secondary anisotropies are probably the most 



40 



Local model 




5 - 
4^ 

3 - 
2- 

1 - 



WMAP (T/T+P) 
Planck (T/T+P) 
CMBPol (T/T+P) 



100' 



50- 
40- 



30- 



200 300 500 700 1000 15002000 3000 



Equilateral model 



20 




WMAP (T/T+P) 
Planck (T/T+P) 
CMBPol (T/T+P) 



200 300 500 700 1000 15002000 3000 



FIG. 11 Fisher matrix forecasts on A/nl, featured for different experiments: WMAP (green, dotted lines), Planck (red, dashed 
lines) and the proposed CMBpol l |Baumann et a/.[[2008| l survey (blue, solid lines). The left panel shows results for the local shape, 
while the right panel refers to the equilateral shape. Thin lines are obtained from temperature data only, and thick lines show the 
improvement in the error bars coming from adding polarization datasets to the analysis for the various experiment. 



important NG sources. Since the main focus on this review is on the primordial bispectrum, we will not describe 
this NG sources in great detail. We will however outline in this section their main effects in order to understand if, 
and how, they could contaminate an estimate of primordial NG. Let us consider a number A^, of sources of a CMB 
bispectrum signal and call B', , ^ the bispectrum produced by the i-th source. Let us also indicate with B^":,^ the 



primordial component of the bispectrum calculated for /nl - 1. For our purposes B^p^^^^' is the signal that we want to 
measure, while the other bispectra are contaminants, that we would like to eliminate. The total bispectrum of the map 
in presence of these contaminants is then: 



1=1 



(nil 30) 



where A, is the amplitude of the /-th bispectrum. If we have a precise prediction of the bispectra generated by the con- 
taminants, we can then think of extending our /nl estimator to a joint-estimator of all the amplitude parameters. The 



optimal cubic /nl estimator defined in (III. 72 1 would then be generalized to the multi-parameter case by minimizing 
the following ;^'^: 



X' 



(/nl,A,) = 2 



(/nl^^p } + 



A W pobs \ 



txhh 



(m.i3i) 



The new errors on /nl in this case can be forecasted as usual by means of a Fisher matrix analysis. The Fisher matrix 
described in the previous paragraph can be generalized straightforwardly to the multi-parameter case. In this case F 
becomes an array whose entries are defined as: 



z 



(III. 132) 



The optimal eiTors on a given amplitude A,- (including /nl) then become, according to the multi-dimensional general- 
ization of the Rao-Cramer bound: 



(nil 33) 
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where the crucial point to notice is that we now first invert the Fisher matrix and then we take the square root of 
the diagonal elements to find the errors. This is the error that is obtained when the full joint-parameter likelihood 
is calculated and then the 1 -dimensional likelihood for a given parameter is obtained by integrating out all the other 
degrees of freedom: a process defined in statistics as marginalization. One can see that the inverse of the Fisher matrix 
defines the covariance matrix of the parameters under study. If the various parameters are completely uncorrelated, 
then the Fisher matrix is diagonal and we would have FJ.^ — {f~^^,,, showing that the parameters can obviously be 
estimated independently and the marginalization process doesn't change the error bars on a given parameter of interest 
(in our case /nl)- If the different parameters are correlated, however, then ofF-diagonal terms appear in the Fisher 
matrix and the error bars after marginalization {i.e. the "real" error bars to quote in the results) are larger than those 
that would have been obtained by naively neglecting contaminants. An obvious but useful observation is that two 
bispectral amplitudes will be strongly correlated when the respective shapes are similar. To make a practical example, 
the bispectrum generated by correlating weak lensing of CMB anisotropics with the Integrated Sachs-Wolfe (ISW) 
eff'ect can be shown to be peaked on squeezed configurations. For this reason the presence of this effect can be a 
significant contaminant for estimates of local non-Gaussianity. 

So far in this section we have described the degradation effects on the error bars if an hypothetical joint-estimator of 
all the CMB bispectrum amplitudes was built, and the amplitudes of contaminants were marginalized over to estimate 
/nl- However a joint-estimation might be difficult, due to factors like the presence of theoretical uncertainties on 
the shapes of contaminant bispectra or possible practical difficulties in finding an efficient implementation of this full 
bispectrum-likelihood estimator {e.g. if the additional secondary bispectra are non-separable). As a result, the practical 
approach so far has been to estimate only /nl using the techniques described in previous sections, and neglect possible 
non-primordial contaminants. In this case the possible effect of contaminants would not show up as a degradation of 
the error bars but in an even worse way, by introducing a bias in the /sjl measurements. Let us see this by assuming 
that the CMB three point function takes contributions both from a primordial NG component and from a contaminant 
bispectrum with amplitude A,. Let us also assume that we can produce a set of NG Monte Carlo simulations of 
CMB maps including both bispectra. We assign a given /nl in input to the primordial component of our simulated 
maps. Finally we estimate the average /nl obtained from the simulations by applying the usual optimal cubic statistic 
described so far. The result of our MC average will be: 

1 d/nl-1 Tiob served 

I f \ - 1 V f|6f3°flf2^3 

where B"'""'"'^'^ is the averaged bispectrum extracted from the map. The /nl term on the r.h.s. of the second line comes 
from the fact that the normalization TV is chosen in such a way as to obtain an unbiased estimator of the primordial 
component. However a second term is present, which accounts for the fact that a contaminant bispectrum, B'^^J'j^ is in 
the map; this term clearly biases the estimator'^. The magnitude of the bias will depend again on how similar the shape 
of the contaminant bispectrum is to the primordial one. If, for example, the contaminant bispectrum is strongly peaked 
on equilateral configurations and suppressed o n squeez ed ones, a local estimator of NG will then not be significantly 



biased by it, since the second term in equation ( 111. 1 34 1 will cancel-out. However, an estimate of f^j^"' will in this case 
be significantly biased. 



Let us remind that by definition an estimator of a parameter A (in our case /nl) is unbiased if (^A^ = A,A being the estimate from data, and A being 
the true parameter of the underlying model. 
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In general we can define the correlation coefficients between two bispectra, labeled / and j, as: 

' . (nil 35) 



The definition of "correlation coefficient" becomes completely transparent if we rewrite the previous formula in terms 
of the Fisher matrix, and keep in mind that F ' define the covariance matrix of the bispectrum amplitudes: 



in. 



(nil 36) 



The correlation coefficient vary by definition from 0, for totally uncorrected shapes, to 1, for identical shapes. The 
more a given contaminant bispectrum is correlated to the primordial bispectrum that we want to measure, the larger 
will be the induced bias. At this point we distinguish between three possibilities. The first is that the contaminant 
bispectrum shape and amplitude are perfectly known. In that case we can compute the expected bias from formula 



(IV.173I and subtract from our estimate. The second possibility is that the shape of the contaminant bispectrum is 
known, but its amplitude is defined with a given uncertainty. In this case we can propagate this uncertainty by quoting 
it in addition to the statistical error bars on /nl obtained in the usual way. The third and worst possibility is that we are 
unaware of the presence of some contaminant effect, or we know nothing about its bispectrum. In this case we might 
obtain a biased estimate of /nl without knowing it and thus eventually misinterpret a spurious NG effect as primordial 
NG. Contaminants are then very dangerous, because if not properly taken into account can lead to spurious claim of 
detection of primordial NG. For this reason, if a positive detection of /nl were to be made at some point for a certain 
model, all possible tests for the presence of contaminant effects should be performed. Moreover, since we cannot 
be absolutely sure that we are considering all possible source of NG contamination, cross-validation of the result 
using other non-bispectrum based estimators will be very important. These other estimators (Minkowski Functionals, 
wavelets, needlets, higher order correlators are just some examples among those considered in the literature) are by 
construction sub-optimal estimators of the primordial component. However, in principle they are expected to produce 
a totally different response to NG contaminants than the primordial bispectrum. A cross-detection of /nl with many 
different statistics would then be much less likely due to some unknown spurious effect. Another way to test the 
primordial origin of an observed NG signal, recently proposed in Munshi & Heavens ( |2009| l, is to modify the optimal 



bispectrum estimator in order to evaluate a function of £ rather than a single amplitude fyi^. The point is that if a 
clear detection of /nl is ach ieved at several o", then the signal is large enough to allow a less radical data compression. 



( 2009 1 have then recently proposed to estimate the "bispectrum related power spectrum", C 



skew 



Munshi & Heavens 
defined as: 

Like in the usual /nl estimator, the optimal S /N weighting is included, and the observed bispectrum from the map 
is correlated to the theoretical shape. However in this case we don't measure the overall amplitude, but rather the 
amplitude for each ^-bin. Note that 

/nl^ iVcf™'; (III. 138) 

by construction of Cf™, the usual /nl estimator is then retrieved by summing the bispectrum-related power spectrum 
over all the i. The general idea is now that the functional dependence of this skew power spectrum on { will show 
significant variation between different sources of NG, allowing a clearer test of the hypothesis that the origin of the 
observed signal is primordial. A number of investigations of WMAP data have already been performed using this 
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statistic in order to look for primordial and secondary sig nals (|Calabrese et al. 2009 Smidt et al. 2009| l, and related 
pseudo-C/ statistics have been developed in Munshi et al. (2009 1. 

In any case, as long as bispectrum estimators are considered, independently of the specific statistic or implemen- 
tation, the best way to deal with NG contaminants is to make sure to list all of them and study their bispectra, or at 
least find ways to assess their potential impact on the final results. In the following paragraphs we will then turn our 
attention to a classification of the most important potential sources of spurious NG, and see how they are treated in 
the primordial bispectrum analysis. Finally, we will consider some effects that interact with the /nl measurement not 
necessarily by directly producing a secondary bispectrum, but rather by changing the normalization of the estimator 
or by increasing the error bars without producing any bias. 



1 . Diffuse foreground emission 



There are three main astrophysical effects producing a galactic microwave emission from our galaxy in the typical 
frequency range of a CMB experiment (Bennett et al. 2003 Dodelson 2003 1: free-free emission from electron- 



ion scattering, synchrotron emission from acceleration of cosmic ray electrons in magnetic fields, and thermal dust 
emission. 

Since these sources produce signals with a peculiar spectral and spatial distribution, multifrequency observations 
allows the separation of them from the primordial component of the CMB signal by suitable component separation 
algorithms. In the resulting "cleaned" map the foreground contribution to the aim is minimized, although obviously it 
can never be completely eliminated. The remaining foreground contamination after cleaning is called the foreground 
residual. Note that the emission from the galactic plane of the CMB map is so strong that a clean separation of the 
primordial CMB component from the foregrounds is impossible. The galactic regions that are too contaminated to 
produce a clean component separation have to be masked out in the analysis. The size of the galactic mask will 
depend on the choice of the foreground flux level above which the pixel is considered too contaminated to be included 
in the analysis. The choice of the cut-off will depend on the specific analysis that one wants to perform on the data. 
Since the primordial NG signal is much smaller than the Gaussian component, more conservative masks (i.e. larger) 
need generally to be used for /nl estimates than those applied to Q estimation. Direct information about the spatial 
distribution of foreground emission in the sky (i.e. free-free, synchrotron or dust) is provided in the form of templates, 
obtained either from the most foreground contaminated channels of the CMB experiment itself, or from external 
astrophysical surveys [e.g. observations of radio-emission, maps of Ho- emission). Templates are affected by several 
sources of uncertainties and errors (see e.g. Bennett et al. 2003') and using them in assessing the possible impact on 
/nl of foreground emission or residuals has both advantages and disadvantages. The safest approach is probably to 
combine internal consistency tests on the data with analysis involving the use of templates. 

The first extensive tests of possible foreground contamination in /nl measurements were performed in Yadav & 
Wandelt| ( |2008 [), where a detection of a primordial local signal at above 99.5% level on WMAP 3-years data was 
claimed. As explained earlier (see caption of table IIIi, further analysis on more recent datasets and/or using more 
optimal estimators have led to an updated f^^- estimate that is about 2-cr away from the origin, i.e. just a "hint" of a 
possible local signal, rather than a detection. However, as long as a detection was claimed in Yadav & Wandelt ( 2008[ l, 
tests to exclude a possible contamination from diffuse foregrounds had to be carried out. In this case the authors relied 
mostly on the "internal consistency test" approach. Their analysis included: 

1. Expanding the original galactic mask in order to see if the estimated value of /nl is stable for different choice 
of the mask. A significantly lower value of /nl for a larger mask might mean that some unmasked noise 
contribution is affecting the measurement with the original mask. 

2. Comparing /nl estimates from foreground reduced maps to estimates from "raw" maps that include a galactic 
mask, but have not gone through a component separation process. If foregrounds have a significant impact on 
/nl then one expects the measurements from raw and reduced maps to differ significantly. 



3. Comparing different frequency channels. If foregrounds significantly contaminate measurements at given fre- 
quencies, then different channels should produce different results. 
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Analysis involving some kind of prior information about foreground emission were carried on both in |Yadav & Wandelt] 
(2008 1 and Smith et al. ( 2009[ ). The two approaches adopted in this case were: 



Producing simulations including both a Gaussian primordial CMB signal and the foreground emission. The 
latter has in this case to be generated according to a model that allows for a good reconstruction of the observed 
templates. The /nl estimator can then be applied to these simulations in order to check if the measured /nl is 
consistent with (as it should be, in absence of significant foreground contamination, since the primordial input 
is Gaussian). 



2. For an optimal estimator including full C pre-filtering ( Smith et al. 



2009 1, adding the foreground templates 
to the noise covariance, by assigning infinite variance to each template T'^^^^^{n). In this way the estimate is 
"blind" to the template amplitudes. This produces a loss of information that in turn determines an increase of 
the variance. The larger the contamination from foreground is, the more the variance increase. For negligible 
contamination, the variance stays the same. In any case, the effect of foregrounds is entirely included in the error 
bars, provided the assumed templates are accurate enough. This method of analysis, called template marginal- 



Rybicki & Pressl ( |1992| ). 



ization, is adopted in Smith et al. ( 2009 1. A complete mathematical derivation of this method is provided in 



In addition to the methods outlined above, there is also the possibility of using the foreground templates for a joint- 

This approach has been recently used 

m 



estimation of /nl and of the templates amplitudes (see equation (III. 131 1). 

Cabella et al. ( 2009| l for a needlet estimator It could be obviously reapplied in the same form to a bispectrum 



estimator. 

In conclusion, all the tests above have been applied to WMAP 3-years and 5-years data releases. No evidence for 
the presence of a significant contamination of the local /nl measurement from diffuse foreground was produced. Other 
shapes of /nl were not considered since the only type of non-Gaussianity that has produced a marginal detection is so 
far the local one. Although diffuse foregrounds and foreground residuals do not seem to contaminate primordial NG 
measurements in WMAP, this is not guaranteed to hold true for Planck, due to its much higher sensitivity. 



2. Unresolved point sources 



Extragalactic point sources are the most important foreground at small angular scales (see Wright et al. 2009 1. 
Sources are identified by searching the maps for bright spots that fit the beam profile, and then masked out. However 
not all the sources can be resolved and eliminated in this way. Unresolved point sources contaminate the map and are 
a source of a NG signal that can potentially interfere with primordial NG measurements. Unclustered extra-galactic 
point sources have a Poisson distribution and their bispectrum is then simply a constant: 



(HI. 139) 



with an amplitude that has to be estimated from the data and depends on the level of contamination from unresolved 
sources. We can now use equation (III.136i to estimate the correlation between primordial shapes and the point 
source bispectrum. For a given choice of the amplitude we can also estimate the expected bias on the /nl estimator 
Simulations of NG maps including the bispectrum from point sources can also be produced and the primordial /nl 
estimator for different shapes applied to them in order to estimate the bias. Finally, since b^,is manifestly separable. 



an estimator of can be built. All of these analysis were performed in Komatsu et al. ( 2009 2003 1 on local and 
equilateral shapes to conclude that point sources do not contaminate significantly the estimate c 

hand they have a larger impact on /^?"'': their induced bias from MC simulations is A/^j"'' - 22 



to the statistical error bar A/^^"" ~ 100. Additional tests were performed in 
possible presence of clustered unresolved point sources. No significant contamination on f^- was found in this case. 
As for the diffuse foreground case, the enhanced /nl sensitivity that Planck can achieve with respect to WMAP might 
increase the impact of these effects. 



Smith et al. 
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3. Secondary anisotropies 



One big advantage of using CMB anisotropies to test primordial NG is that they are small and can then be treated 
in the linear regime. The CMB temperature fluctuation field is thus linked to the primordial potential through a 
linear convolution with radiation transfer functions, as we saw earlier. At this level, the Gaussianity of the primordial 
potential is conserved in the CMB temperature fluctuation field. If, however, we work at second order in perturbation 
theory, the initial conditions are propagated non-linearly into the observed CMB anisotropies, and the resulting CMB 
fluctuations are mildly non-Gaussian even starting from a Gaussian primordial curvature field. Second order effects 
are clearly very small. However they may well be of the same order of magnitude as primordial NG, since the NG 
component of the primordial potential is O (/nl ('^l^'^)))- conclusion, secondary anisotropies are a potential source 
of CMB NG, at a level that could in principle contaminate estimates of primordial non-Gaussianity. To fully account 



for these effects, it is necessary to obtain a relation analogous to equation (III.30i, but to second order in perturbation 
theory. Radiation transfer functions are obtained at first order by solving the linearized system of Boltzmann-Einstein 
equations (see e.g. Dodelson, 2003, Ma & Bertschinger 1995. ). The same equations will then have to be expanded 
and numerically integrated at second order in this case. Having obtained second order transfer functions, the full 
angular bispectrum of secondary anisotropies can be calculated and correlated to the primordial one in order to check 
for the presence of contaminant effects. A full second order Boltzmann code is actually not yet available, although 
much progress has been made over the past few years in that direction. The full system of second order Einstein- 
Boltzmann equations has been derive d in ([Bartolo et al. 2006a|b Pitrou 2009 Pitrou et al. 2008 ) and partially 



integrated numerically in ( |Nitta et oT] |2009| l including only the source terms that can be written as product of first 
order perturbations. These terms have been shown to produce a totally negligible NG contamination. No numerical 
evaluation is available to date of the "genuine" second order source terms. Although a full solution of the relevant 
equations hasn't been obtained yet, a significant number of secondary effects are known, and have been modeled 
for some time. Among these there are for example weak lensing, Sunyaev-Zeldovich (SZ) effect, Rees-Sciama (RS) 
effect, and so on. Therefore, a natural approach that was adopted in the literature consisted in studying the bispectra 
arising from these well-known effects and from their correlations {e.g. ISW-lensing correlation, SZ-lensing correlation 
and so on). It goes beyond the purpose of this review to discuss in detail these results and their implications. Let us 
just mention them briefly. A fisher matrix analysis in |Serra & Cooray ( 2008 ) showed that the combination of bispectra 
arising from ISW-lensing, SZ-lensing and unresolved point sources produced a negligible contamination at the angular 
resolution and sensitivity of WMAP, but a significant one for an experiment with the characteristics of Planck. It was 
in particular shown that estimates of local NG would be biased, especially by ISW-lensing correlation, with /j^j"* ~ 10 



for local NG. A similar result on ISW-lensing was obtained in another Fisher matrix analysis by |Smith^ Zaldarriaga 
( 2006) 1, and a similar level of contamination was found in Mangilh & Verde (20091 for the analogous RS-lensing 
bispectrum. A bispectrum estimator of local and equilateral NG was applied to simulated lensed primordial NG CMB 



maps by Hanson et al. ( 2009 1, and three main effects were studied: a possible bias induced by neglecting the lensing 
of primordial bispectrum in the normalization and weights of the estimator, an increase of the variance due to lensing- 
produced higher order correlators, and ISW-lensing bias. The only significant effect turned out to be the ISW-lensing 
bias on /^"l , at a level confirming Fisher matrix predictions. Note that this bias, being well-known and expected, can 
be simply calculated and subtracted from future Planck estimates. The reason why the coupling between lensing and 
ISW tend to bias the local estimate can be understood physically: large scale potential fluctuations source the ISW 
effect and produce a lensing signal on small scales, generating a NG signal on squeezed triangles. Although both the 
primordial local bispectrum and the ISW-lensing bispectrum are peaked on squeezed triangles, the presence of acoustic 
oscillations in the primordial configurations reduce the overall correlation between the two shapes, thus making the 
final bias significant, but not too large. In order to conclude our brief survey of studies of secondary bispectra, let us 



mention the work done in Babich & PierpaoU ( 2008 1, where point source density modulation bispectra induced by 
lensing magnification and selection effects, as well as SZ modulation from lensing magnification were studied. The 
conclusion was again that these effect are negligible for WMAP but close to the sensitivity level of Planck. Despite the 
great attention received so far in the literature, much work still has to be done in the area of assessing NG contamination 
from secondary sources. It is clear that a complete and accurate description of secondary bispectra will be crucial for 
analysis of the future Planck data set. 
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4. Non-Gaussian noise 

Systematics are another potential cause of contamination beyond astrophysical and cosmological sources. The 
noise in the experiment is generally well described as Gaussian. However possible non-Gaussian properties have to 
be tested in our context. This was done in |Yadav & Wandelt ( 2008 1 by taking differences of yearly WMAP data in 
order to create jackknife realizations of WMAP noise maps for different detectors, including instrument systematics. 
The estimator can then be applied to these realizations in order to check that a negligible /nl is measured. This was 
the result obtained on the WMAP 3-year data-set. 



5. Other effects 



In this section we quickly summarize other effects that could interfere with estimates of primordial non-Gaussianity, 
but did not fit the classification above in the sense that they do not correspond to NG effects contaminating the CMB 
sky or the instrument noise. 

One of this effects is 1 // noise, expected to affect especially the low frequency channels of Planck. The 1 // noise 



component is generally removed from the map using "destriping" algorithms (see e.g. Efstathiou 2005 Maino et al. 



20021. The unsubtracted "destriping residuals" form a Gaussian correlated random field in pixel space. Their non- 
trivial covariance matrix should in principle be included in the inverse covariance pre-filtering of the optimal estimator 
If not included in the pre-filtering, this effect could in principle enhance the estimator error bars (although it cannot 
generate any bias, since it is Gaussian). Unfortunately, a full numerical evaluation of this covariance matrix is quite 
challenging. Donzelli et aZ.|(|2009| applied the estimator in its pseudo-optimal implementation to maps of Gaussian 
CMB signal + noise, accounting only for anisotropic noise in the linear term, but including destriping residuals in 
the noise model adopted for the simulations. The final result shows that the error bars do not increase when 1// 
noise effects are included in the simulations, even though they are neglected in the covariance matrix appearing in the 
estimator. 

Another effect to take into account for Planck is that of an asymmetric beam. The beam in the estimator normal- 
ization term is approximated as a circular beam. However Planck optical simulations (e.g. Sandri et al. 2004 1 show 



that in reality we have to deal with elliptic beams, characterized by a non-trivial azimuthal dependence. If the circular 
beam approximation in the normalization of the estimate is not accurate enough, a bias could be introduced. Moreover 
the anisotropy of the beam could cause an increase of the variance if neglected in the inverse covariance pre-filtering. 



Again, these effects were found to be negligible in tests on realistic simulations performed by Donzelli et al. (2009 1. 

Finally, the estimate of /nl is done assuming a given cosmological model, i.e. by fixing all the other cosmological 
parameters to their best-fit value obtained from a likelihood analysis of the angular power spectrum. Since they are 
themselves the product of a statistical estimation process, these values obviously present uncerta inties that should be 
propagated into the final /nl -error bars'"^. This calculation was done in Liguori & Riotto (2008 1, where it was found 



that the propagated error is /ml dependent and it can become important only if a large /nl will be detected in the data 
at some point. 



J. Generation of simulated non-Gaussian CMB maps 

In this section we will describe algorithms for the generation of non-Gaussian CMB maps with a given bispectrum. 
There are three main reasons why primordial NG simulations of the CMB are useful in the context of bispectrum 
estimation of /jl: 



In particular, since we are not doing a joint-likelihood estimation of all the parameters and marginalizing to get /nl (that would be the optimal but 
time consuming approach), the effect of uncertainties in the parameters propagate onto the /nl measure as a fora.v. This bias has to be evaluated 
and quoted in addition to the usual statistical /nl error bar 
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FIG. 12 Top panels: a Gaussian realization of the CMB sky (fe /f) and a non-Gaussia n local CMB map (right), obtained by adding 
to the Gaussian one a NG component with f^^- = 3000. FromlLiguori ef a/. J2003|. L ower panels: a Gaussian CMB map (left) 
and a non-Gaussian DBI map (right) with f^^' = 4000. From (Fergusson et at. 3009 l The maps in the upper panel have been 
obtained using the local algorithm described in section [ni. J. 1 1 The maps in the lower panel have been produced with the bispectrum 
algorithm of section [lII.J.2| after having separated the primordial DBI shape using the eigenmode expansion defined in ( |II.26^ . 



1. To test the unbiasedness of the /nl bispectram estimator (by checking that the Monte Carlo average of the 
recovered /nl reproduces the /nl set in input). 

2. To study how the expected primordial NG signal imprinted in the CMB is modified by the presence of other 
effects, like those considered in section For example, weak lensing of primordial NG might in principle 
change the observed bispectrum and affect the estimates. This can be studied again by testing the estimator on 
NG lensed simulations, as it was done in |Hanson et a/.| ( |2009| l 

3. For local NG, to obtain the error bars of the /nl estimator if a large /nl is detected at several cr (see section 



III.E.4 1. We have previously seen that for a several sigma detection of local NG the bispectrum variance is /nl- 



dependent. The Monte Carlo average ( III. 1 17 1 thus has to be evaluated on NG simulations with the measured 
/nl in input. 

Unless we are in the situation described at point 3. of the list above all we need to produce is then maps with given 
power spectrum and bispectrum, since higher order correlators can be neglected. In the large local /nl case higher 
order correlators are instead important and have to be included. Fortunately the local case is the only one for which 
we have a full expression of the primordial potential <i>(x) that allows us to produce exact simulations. 

We will divide this section in two parts. In the first we will describe exact simulation algorithms of local NG, while 
in the second we will describe methods to generate maps with given power spectrum and bispectrum, starting from an 
arbitrary primordial shape. 
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1 . Algorithms for local non-Gaussianlty 

First of all, let us recall that the CMB multipoles aim are related to the primordial gravitational potential O through 
the well known formula: 



where Af(fe) are the radiation transfer function and the potential is written in Fourier space. We already met this 



«fm= r ^H^)Yt„r{k)ht{k) , (HI. 140) 



formula when we calculated the relation between the primordial and CMB bispectrum in section III. A We also recall 



that the local non-Gaussian primordial potential takes a very simple expression in real space, where: 

<l)(x) = O^x) + /nl [<l>i(x) - (<l>i(x))] . (ni.l41) 

In the previous expression Ol is a Gaussian random field, characterized by a primordial power spectrum f o(^) = 
Afc" in the following we will refer to <1>l(x) as the Gaussian part of the primordial potential. The remaining non- 
Gaussian part of the potential is simply the square of the Gaussian part point by point (modulo a constant term, 
necessary to enforce the condition ('l'(x)) - 0; however it is clear that this term only affects the CMB monopole). It is 



then convenient to work directly in real space and recast formula ( III. 140 1 in the following form 

/ 



atm^ \ d\^(r)Yt,n{r)ae{r), (III.142) 



where a{{r) = J dkk^j((kr)Af(k), also used in (III.31 1, is the real space counterpart of the radiation transfer functions 
A({k), jt{kr) is a spherical Bessel function, and r is a look-back conformal distance. This formula suggests to structure 
an algorithm for the generation of local CMB NG maps in the following steps 

1. Generate the Gaussian part Q>i of the potential in a box whose side is the present cosmic horizon. 

2. Square the Gaussian part point by point to get the non-Gaussian part. 

3. Expand in spherical harmonics the Gaussian and non-Gaussian parts of the potential for different values of the 
radial coordinate r in the simulation box. 

4. Convolve the spherical harmonic expansions of <1>l and Onl with the radiation transfer function Af(r) in order 
to obtain the Gaussian and non-Gaussian part of the multipoles of the final NG CMB simulation. For a given 
choice of the non-Gaussian parameter /nl a CMB map is then obtained simply through the linear combination 
at„, - a^^^ + /nla^ (the superscripts L and NL always indicating Gaussian and non-Gaussian respectively). 

The most difficult and time consuming part in this process is actually the generation of the Gaussian part of the 
potential O. One possibility is to generate the Gaussian part of the potential in a cubic box in Fourier space (where 
different modes are uncorrected and have variance given by the primordial power spectrum P<^(k), then apply a Fast 
Fourier Transform (FFT) algorithm to go to real space. Cartesian coordinates are then transformed into spherical 
coordinates by means of an interpolation algorithm in order to transform <l)/.(x) into <^L{f, «)■ Finally, the Gaussian 
potential in spherical coordinates is squared point by point to get the NG part and the spherical harmonic expansion 
and radiation transfer function convolution at point 4 of the list above are performed in order to obtain the multipoles 



of the final CMB map. The aforementioned algorithm was implemented in Komatsu et al. ( 2003 1 to generate NG local 
CMB maps at the resolution of the Planck satellite. 

The difficulty with this approach arises from the fact that we are working in a box of the size of the present cosmic 
horizon, (about 15 Gpc in conformal time), but at the same time a cell in this box must have a side no bigger than 20 
Mpc in order to resolve the last scattering surface, where most of the CMB signal is generated. A more convenient and 



accurate way to produce the local NG flf,„ was found in (Liguori et al. . 2006 ; Liguori et al. 2003]l: the idea is to work 
directly in spherical coordinates, use a non uniform discretization of the simulation box (since no sample points are 
needed in a large region of the box where photons are just free streaming, while many sample points are needed at last 
scattering, as we just pointed out above) and generate the multipoles of the expansion of <1>l(x) through the following 
two step approach: 
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1. Generated uncorrelated radial multipoles n(m(r), gaussianly distributed and characterized by the following spec- 
trum: 

(«.„„(r,)4„,(r2)> = ^"^y'^ 4]C] ; (ni.i43) 

where 6d is the Dirac delta function. 

2. Filter the multipoles nf,„ with suitable functions in order to produce a Gaussian random field with the properties 
of the multipole expansion of the primordial Gaussian potential . It can be shown that the expression of the 
filter functions is: 

Weir, ri)^^Jdkk^ je(kr)ji(kri ) , (III. 144) 

where Po is the primordial curvature power spectrum, and the filtering operation takes the form 

'^'LW^ jdrirlne„,{ri)We{r,ri) . (III.145) 

In the last expression <t^^j^(r) are the desired quantities, i.e. the multipoles of the expansion of the Gaussian part 
of the primordial potential for a given r. 

This algorithm, recently improved in |Elsner & WandeTt ( 2009 1, was used to produce NG local maps at the resolution 



of WMAP and Planck in temperature and polarization. An example of its results is shown in the upper panels of fig. 12 
Fig. 13 shows 1-point PDFs of temperature anisotropics for different values of % extracted from these simulations. 



2. Algorithms for arbitrary bispectra 



In the limit of weak non-Gaussianity, an algorithm to produce non-Gaussian CMB simulations with a given power 



spectrum and bispectmm for separable primordial shapes was described in Smith & Zaldarriaga ( 2006) . In this algo- 
rithm the non-Gaussian components of the CMB multipoles are obtained using the following formula: 



■'em 



e,m, ^ 



m2 



t3 

m3 



(III. 146) 



where a^^^ is the Gaussian part of the CMB multipoles, generated using the angular power spectrum Cf, while Be^e^jj 
is the given bispectrum of the theoretical model for which simulations are required. Note that alternative algorithms to 
generate CMB maps with given bispectrum have been proposed in the literature (Contaldi & Magueijo , 200 1 ; Rocha] 
et al. 2005 1 ) but they are less general that the one introduced by equation ( III. 146 1. Although ( III. 146 1 is completely 
general, as before its numerical evaluation is only computationally affordable for bispectra that can be written in 
separable form. W e have emphasized already that separability results in a reduction of the computational co st of the 
estimator (III. 100 1 from 0({f„gj to <9(^^,ax) operations; the same argument applies here, and allows to rewrite (111.146 1 



into an equivalent form in pixel space. Starting from formula (III. 46 1, and substituting it in (III. 146 1, we find: 



dQ.f, (2Xe(r)MY{r, n)Mz{r, n) + 2Ye{r)Mx(r, n)Mz(r, n) + 2Ze(r)Mx(r, njMyir, n)) . (ni.147) 



As aheady discussed in the /NL-estimator section, the limitation dictated by separability is clearly overcome by using 
the eigenfunction representations for the bispectrum (11.26 1 and (III.6O1 introduced in Fergusson et al. (2009 1. As 



usual, the basic idea is to start by expanding an arbitrary bispectrum shape S (either primordial or in the CMB) using 
a separable polynomial decomposition until a good level of convergence is achieved and then to substitute the mode 



decomposition into (III. 146 1 to get a linear combination of numerically tractable terms written in the form ( III. 147 1. 
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Probability density function of tempera ture pi xel from local primordial non-Gaussian CMB maps, obtained with the 



FIG. 13 

"exact" simulation algorithm described in section |III.J.l 



Different panels show the result for different values of f^^ , in order to 
give an idea of the order of magnitude of the signal that one wants to detect. For /nl < 1000 the non-Gaussianity is too small to be 
seen this this plots. Note that WMAP constrain f^[- to be < 100. 



Using the separable mode coefficients ap, s for the reduced bispectrum ( III. 56 1 and the filtered map expressions Mp(r, h) 
(III.llOi as the starting point, we find that the expression (III.147i generalises to 



" A Z "i"--' j dxxWp(x) j d^h YTi^) M^ir, n) (r, n) , 



(HI. 148) 



where the M^(r, n) are found by summing using a set of Gaussian aS^'s convolved with the q' ^s (refer to eqn (III. 108 1) 



M' 



Im 



Ci 



(HI. 149) 



Here, the accuracy of convergence with the Up^s is parametrized in terms of the correlation C(S,Sn) between the 
original non-separable shape and the eigenmode expansion, as defined previously (11.19 1. Note that this convergence 
can also be checked more accurately using the full Fisher matrix coiTelation on the CMB bispectra C{btit^e3,b'^ { ( ), 
described in sections UlI.HI and I lII.II 

In addition to the bispectrum separability req uirement , there is an important further caveat which can prevent the 
straightforward implementation of the algorithm (III. 146 1. By construction, terms 0(f^^ and higher are not explicitly 
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controlled. Following the discussion in (Hanson et al. 2009) we can write the connected N-point functions as: 

= [C^, +/NLCrl (HI. 150) 



(at, 



(111.151) 

(ni.i52) 



Thus the condition that the map has the pow er spectrum C/ specified in input will only be satisfied if the power spectrum 



of the non-Gaussian component in ( 111. 150 1 remains small. Since this method does not control 0(fc,, ) terms, one has to 



ascertain that spuriously large C™ contributions do not affect the overall power spectrum significantly. It turns out that 
this effect plagues current map simulations if the standard separable expressions for the local and equilateral bispectra 



are directly substituted into (III. 146 1. However a slight modification of equation (111.146 1, described in Hanson et al. 
([2009 ) and Fergusson et al.' (2009) allows us to overcome this problem at no computational cost. Moreover, it was 
shown by Fergusson et a/., ( ,2009) that maps obtained from the eigenmode expansions (11.26 1 and ( 111.60 1 are stable 
independently of the shape under study, thus making this map-making generating algorithm robust and fully general. 
Examples of DBl NG maps produced by combining the eigenmode expansion method with the map making algorithm 
described in this section are shown in the lower panels of fig. [12] 
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IV. LARGE-SCALE STRUCTURE 

In the standard scenario, early perturbations produced during inflation are responsible for the common origin of the 
CMB temperature fluctuations and the large-scale matter and galaxy distributions in the Universe, i.e. the large-scale 
structure. The Cosmic Microwave Background provides a remarkable example of a Gaussian random field in nature. 
Information on cosmological parameters is in fact derived from measurements of its power spectrum, the C/'s, while 
bispectrum measurements from WMAP data remain consistent with zero. The distribution of matter, as we can infer 
today from shear or galaxy observations, unlike the CMB, can be described as a highly non-Gaussian random field, 
even for Gaussian initial conditions. 

The matter overdensity 5(x) is defined in terms of the matter density p(x) and its mean value p by 

^(x) = P^^'^'P ^ (IV. 153) 

P 

with zero mean by construction. Here, at late times, the limiting value 5 - - I 'm voids, accounting for a large fraction 
of the volume of the Universe, while achieving values 5 » 1 in collapsed objects such as dark matter halos. Its 
probability distribution function is therefore expected, at low redshift, to depart strongly from a Gaussian distribution 
centred at 5 = 0, even though it could be well approximated by it at decoupling, when perturbations around 5 - Q 
were of the order of 5 ~ 10"^. Such non-Gaussianity is the result of the nonlinear evolution of structures subject to 
gravitational instability. 

In addition, nonlinearities in the bias relation between the galaxy and matter distributions constitute a second source 
of non-Gaussianity in the large-scale structure mapped out by redshift surveys. Non-Gaussian initial conditions would 
therefore provide a third component in the non-Gaussianity of the galaxy distribution. The question regarding the 
detection of effects due to primordial non-Gaussianity, is therefore strictly related to our ability to distinguish between 
these difi'erent contributions and, ultimately, it will depend on the robustness of our theoretical predictions in the linear 
and mildly nonlinear regime. From this respect, cosmological Perturbation Theory (FT), and its more recent develop- 
ments, is very important for providing the tools to study the evolution of non-Gaussianities and how to differentiate 
their origin. 

Considering only the matter distribution, the leading order prediction in standard PT for the matter bispectrum at 
large scales is given by the sum of a primordial component and a component due to gravitational instability, which 
is present also for Gaussian initial conditions. Until fairly recently it was assumed that this picture could be easily 
extended to the galaxy distribution, with the galaxy bispectrum receiving an additional contribution due to nonlinear 
bias. Following the historical development of the subject, in section [lV.A| we will discuss early work on higher-order 
moments of the matter and galaxy distribution, starting with the skewness. Here, most of the theoretical results on 



higher-order correlation functions are developed. We will then consider in IV.B the matter bispectrum and its de- 
scription in Eulerian perturbation theory, with specific attention given to effects at large scales due to a primordial 
component, as well as at small-scales, nonlinear corrections in presence of non-Gaussian initial conditions. In sec- 
tion IV.C we will deal with the galaxy bispectrum. We will first introduce the simple model based on local bias and 
discuss problems related to bispectrum measurements in redshift surveys with specific attention given to the detec- 
tion of primordial non-Gaussianity. We will see how early results indicated that the galaxy bispectrum could be used 
as a tool to constrain non-Gaussian initial conditions which is, in principle, competitive with the CMB, illustrating 
this with actual results from current data-sets. We will then consider the outcome of recent N-body simulations with 
non-Gaussian initial conditions showing that the simple prediction for the galaxy bispectrum assumed in most of the 
previous literature on the subject fails to describe not only the measured halo bispectrum, but even the halo power 
spectrum, even at large scalesl We now know that correlators of biased populations such as galaxies and dark matter 
halos, receive large corrections, at large scales, from local primordial non-Gaussianity. These results opened up new 
and promising opportunities for detection in future large-scale structure observations. Although, in our view, a proper 
understanding of these efifects remains to be adequately developed at the time of writing, particularly with respect to 
higher-order galaxy correlation functions, we will describe the different descriptions proposed so far in the literature 
and the prospects for detection of primordial non-Gaussianity in measurements of the galaxy bispectrum. 

From an historical perspective, non-Gaussian initial conditions have been studied for quite a long time. For instance, 
early works on the clustering of density peaks and rare objects can be found in |Grinstein & Wi^ ( |1986l l; |Lucchin| 
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& Matarrese ( 1988 1; Matarrese et al. (19861, while early N-body simulations with non-Gaussian initial conditions 



go back to the early eighties (Coles ef a/. 1993 Messina ef a/. 1990 Moscardini ef a/. 1991 Weinberg & Cole 



1992 White 1999 1. In the early days, a large variety of non-Gaussian models, often defined in terms of a nonlinear 



transformation of a Gaussian field were considered. In some cases, a large non-Gaussian component were studied 
because, on the one hand, they could be used falsify some models and, on the other, as a way to reconcile contradictory 
observational results with theoretical frameworks. In this review, however, we will consider only models predicting 
small departures from Gaussian initial conditions which are consistent with CMB observations. 

Finally, while we focus in this review on direct bispectrum measurements, it should be stressed that the effects 
of primordial non-Gaussianity on large-scale structure are not limited to corrections to its higher-order correlation 
functions. Aside from the recent results on the galaxy power spectrum mentioned above, significant departures from 
Gaussian initial conditions are expected to have important effects on the halo mass function and therefore on the 
observed cluster number density. See section 2.1 in Sefusatti et al. ( 2007[) for an brief overview of previous work 
anc^Afshord i & Toll ey|(p008); [Dala i et a/.'('2008^;'Desjacques et a/.'C2009');'FedeH et a/.'('2009^;'Grossi et a/.'('2009l); 



Lam & Sh eth (2009);" lLoVerde etal. (2008 ); Maggiore & Riotto (2009 ); Oguri (2009); Pillepich et al. (2008); Valageas] 
(|2009) for recent theoretical and N-body results. In addition, the corresponding effect on the abundance of voids has 
been studied by |Kamionkowski et al.\ ( |2009) l, while the possibility of constraining primordial non-Gaussianity from 
measurements of Minkowski Functionals in large-scale structure has been explored by Hikage et al. (2008 2006[ l. 



Further effects on the intergalactic medium and reionization (Crociani et al. 2009 Viel et al. '2009) or on future 21cm 



observations ( Cooray 2006 Pillepich et al. 2007| have also been investigated. We refer the reader to other reviews in 
this issue for a more complete discussion of these alternative approaches. 



A. The skewness 



Since the first large-scale observations did not allow an accurate determination of individual bispectrum or trispec- 
trum configurations, most of the attention in the early literature focused on the moments of the galaxy distribution, and, 
in the first place on the third- and fourth-order moments, i.e. the skewness and kurtosis, respectively. The "normalized" 
moment of order p can be defined in terms of the smoothed density field 6r{x) as 



<5^(X))W2 ' 



(IV. 154) 



For Gaussian initial conditions, a perturbative treatment of the equations of gravitational instability predicts at leading 
order ( |PeebIesl[T980l l 

34 



(IV. 155) 



with cr^ - {S]f), computed in Unear theory. When non-Gaussian initial conditions are present, one expects an extra 
contribution to the skewness, typically with a different relation with ctr, whose value depends on the non-Gaussian 
model. Comparisons between the second- and third-order moments, S's r and ctr, (as well as higher-order moments 
such as the kurtosis) measured in redshift surveys have been early recognized as a tool to test the Gaussianity of 



primordial perturbations 


(Bouchet et al. 


1992; Coles & Frenk 


1991; Coles et al. 


1993 


|Fry & Scherrer 


1994 


Juszkiewicz & Bouchet, 


1992, .Juszkiewicz et aZ.,,1993,.Lahav et a/., ,1993, ,Lucchin et al.\ 


19941 |Luo & Schi-amm 


1993 1. These works recognized as well the importance of rehable predictions in the nonlinear regime and of a proper 



modeling of the effects of galaxy bias. In this respect, |Fry & Scherrer ' (1994) proposed a more quantitative prediction 
for the contribution to the galaxy skewness due to galaxy bias based perturbation theory and on the local bias expansion 



of (Fry & Gaztanaga 1993 1. They derived, for the skewness of the galaxy distribution, an expression of the form 

34 



6/72 

7 bi 



(IV. 156) 



where we assumed non-Gaussian initial conditions described by a non-vanishing initial skewness s^"^ (but vanishing 
higher-order moments) and where bi and b2 represent constant bias parameters typical of the galaxy population (which 
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FIG. 14 Left panel: measurements of the skewness of the matter distribution in N-body simulations as a function of the smoothing 
scale R for the non-Gaussian texture model {filled symbols) and for Gaussian initial conditions {open circles). The time evolution 
is parametrized by the value of erg, with the non-Gaussian results shown at erg = 0.4 {triangles) and erg = 1 {squares). Lines show 
different theoretical predictions. Right panel: measurements of the third, fourth and fifth-order moments of the galaxy distribution in 
the APS Galaxy Survey, compared with the simulation results with non-Gaussian initial conditions with dilferent bias assumptions. 
From |Gaztanaga & Mahonen| ( [1996| > (see the reference for further details). 



we will discuss explicitly in section IV.C i. This relatively simple expression describe the skewness measured in galaxy 
surveys, as the sum of three components corresponding to three sources of non-Gaussianity for the galaxy distribution: 
one primordial, one due to gravitational instability and the last to nonlinear bias. Further studies in perturbation theory 
can be found in ( jChodorowski & Bouchet| [1996; Dur rer et al. | [2000 ) while an alternative derivation of the smoothed 
moments of the density field based on the spherical collapse model has been studied in ( [Gaztanaga & Fosalba||1998| . 
The skewness predicted by texture models has been studied in simulations as a function of the smoothing scale R 
by 'Gaztanaga & Mahonen' (T996\ and compared to measurements of the same quantities in the APS Galaxy Survey 
(Gaztanaga, ^ 1 994^ ), see fig. 14 The differences between the 53.R in the non-Gaussian texture model with respect to the 
Gaussian case provides a qualitative example of the typical effects that we expect for non-Gaussian initial conditions 
as a function of the smoothing scale R and redshift. 

The measured skewness, as higher order moments, coiTesponds to a single number Despite the possibility to study 
its peculiar dependence on the smoothing scale R, it is nevertheless difficult to separate the different components, 
particularly with respect to bias effects. However, this possibility is offered in principle by direct measurements of the 
galaxy bispectrum, relying on its dependence on the shape of triangular configurations. In the next sections we will 
discuss in details first the bispectrum of the matter distribution then the bispectrum of the galaxy distribution, a direct 
observable in redshift surveys. 



B. The matter bispectrum 



In this review we will focus on the predictions for coiTelation functions in Fourier space from Euleiian Perturba- 
tion Theory (PT). This approach solves perturb atively the equations for the matter density and velocity field evolution 
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governed by gravitational instability. These are the continuity equation, the Euler equation and Poisson equation re- 
lating the matter density and the gravitational potential. In the PT framework, the relation between the results and 
the initial conditions, given in terms of the initial correlators of the density field is particularly transparent. Moreover, 
recent works have significantly extended, as we will discuss later, the predicting power of this specific tool. Different 
approaches are also available: see, for instance, Scoccimarro (2000a I for a comparison between bispectrum measure- 
ments in N-body simulations and predictions in Lagrangian Perturbation Theory. We refer the reader to Bernardeau 



et al. ( 2002 1 for a comprehensive review of cosmological perturbation theory of the large-scale structure. 



1 . Leading-order results in Perturbation Theory 

As mentioned before, we consider specifically models where non-Gaussian initial conditions are completely given 
in terms of the correlators of the curvature perturbations at early times, and the mechanism responsible for the extra 
non-Gaussian properties of the density field is not active during the subsequent evolution of matter perturbations, 
governed only by gravitational instability. In PT, the so lution for the evolved matter density contrast is expressed as a 



series of corrections to the linear solution (5*'' (Fry 1984j ) 



5k = 4" +4" + ■■■■' (IV. 157) 

where each term can be written formally as'^ 

5^") = J d'qi . . . cPqnFMx , . . . , q„) < . . . , (IV. 158) 

with F„{{{\,. . . , q„) representing the symmetrized n-order kernel in PT. The initial conditions in the Gaussian case are 
completely specified by the linear power spectrum Po(k), with (^j; '5^',^) = SD(^i2)Po(ki), where we adopt the notation 
kjj = kj + kj. Non-Gaussian initial conditions are described, in the first place, by a non-zero expression for the three- 
point function of the linear solution, that is (^ti'^^V'^kj')- In turn, the initial matter correlators, i.e. the correlators of the 
linear solution (5^'\ are given in terms of the correlators of the curvature perturbations, as 

<5k, ■ ■ ■ K) = M{kuz) ■ ■ ■ M(k„,z) <(Dk, ■ ■ ■ 1>k„), (IV.159) 

where we introduce the function 

with T{k) being the matter transfer function and D{z) the growth factor, expressing Poisson's equation in Fourier space 
as 

dviz)^ M{k,z)^v.. (IV. 161) 

Notice that we denote with the primordial curvature perturbations, i.e. evaluated during the matter dominated era, 
not their value linearly extrapolated at present time'^. The linear, i.e. initial, power spectrum is given by 

Po(k) = M^{k, z)P^{k) , (IV.162) 



" From now on we will adopt a different convention for the Fourier transform with respect to the one used for the formulae in previous section. The 
present convention is more common in the large-scale structure literature and conforms with the one adopted in the classical review [Bemardeau] 
I gfa/.|j2002) . 

" This choice, not unique in the literature, is particularly convenient since curvature perturbations are constant during matter domination. Also, it 
conforms to the definition of /nl in terms of <D assumed in the CMB literature on observational constraints and specifically in |Komatsu & Spergel| 

pool) . 
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while the initial bispectram and trispectrum are 

BQ{ki,k2,h) = M{h)M{k2)Mih) B^{kuk2,h) , (IV.163) 
To(ki,k2,k2,k4) = M{k{)M{k2)M{k^)Mik^)T^{ki,k2,h,h)- (IV.164) 

Notice that given these simple relations between curvature and primordial matter correlators, issues such as the prop- 



erty of separability discussed in section II. B for the CMB bispectrum are not present in the case of three-dimensional, 
large-scale structure observables. 

The nonlinear power spectrum is obtained perturbatively from the expansion 

<5k,5k,) = <<'/4?) + + <4?C> + <^k?C> + ■ ■ ■ ' ^^^-165) 

where the term (<5^'*5[;',*) corresponds to the linear solution, while the other terms represent, in analogy with 

perturbation theory in quantum field theory, one- and higher loop corrections as they involve integrations over internal 
momenta. In particular, the term (<5^'^<5^?) vanishes for Gaussian initial conditions as it depends on the initial bispec- 



trum Bo (see Taruya et al. 2008jfor an analysis of nonlinear corrections to the matter power spectrum due to primordial 



non-Gaussianity). 

In a similar fashion, nonlinear corrections in ( |IV.157[ ) provide a perturbative expansion for the matter bispectrum, 

i6,,6,M,) . i^^^) . <C<'<') . <C<^^'C) + iS^yX^ - (I^^166) 

In this case, the leading order contributions are given by the tree-level terms (^^''^'.'^k?) <^k?4?^k?>' wi* the first 
being the initial component and the second corresponding to a contribution to the matter bispectrum due to gravity 
alone, of the form 

B'^''(kuk2, h) = 2F2(ki, k2)Po(^i)i^o(^2) + 2 perm. . (IV.167) 

Notice that this contribution is present even for Gaussian initial conditions as it depends only on the initial power 
spectrum Pq and describe the emergence of non-Gaussianity due to gravitational instability. The leading order, tree- 
level expression of the matter bispectrum with non-Gaussian initial conditions is therefore given in terms of the sum 

B"-''{kuk2M) = B(,(kuk2M\z) + B'[r\kuk2M\z) . (IV.168) 



This expression corresponds to the first two terms on the r.h.s. of (IV. 156 1 for the skewness, which can be obtained 



from (IV.168 1 by integration. 

The possibility of distinguishing the primordial component Bq from the gravity-induced one Be relies on their 
specific and distinct dependence on scale, on the triangular configuration shape and on redshift. For a primordial 
non-Gaussianity described by a curvature bispectrum obeying the hierarchical scaling Z?o ~ typical of weakly 
non-Gaussian models such as the local and equilateral ones, the different redshift and scale dependence of the two 
contributions is evident in their ratio for equilateral triangles {ki - k2 = kj, - k), given by'^ 

Bojk, k, k; z) ^ 1 /nl k-M /nl aV.169) 
B'[.'%k, k, k; z) 4 Mik; z) ~ k^D(z) ' 

We therefore expect, for a wide range of non-Gaussian models, the initial contribution Bq to be larger at large scales and 



at high redshift. The upper left panel of fig. 15 shows the two contributions and their sum for equilateral configurations 



B{k, k, k) as a function of k. The other panels show the effect of the primordial component for different non-Gaussian 



The first equality is in fact identical for local, equilateral and orthogonal non-Gaussianity. simpl y' by definition of the equilateral bispectrum, 
Eq. ??, introduced in I Babich et al. 2004 1 and of the orthogonal bispectrum, Eq. ??, introduced in I Senatore et al. 2009), where f'^ and f^°^' 
are precisely the amplitudes that provide the same value for the curvature bispectrum as the local model tor equilateral configurations. 
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FIG. 15 Effect of the primordial component for different non-Gaussian models on the equilateral configurations of matter bispec- 
trum, B{k,k,k), at redshift z = 1, as a function of scale, at tree-level in PT. In the upper left panel the continuous line shows the 
initial component B(,{dotted line), the gravity-induced component, 5^'^' (dashed line) and their sum (continuous line). For equilat- 
eral configurations the initial component coincides for the local, equilateral and orthogonal models while it vanishes in the folded 
model. In the other panels, continuous lines show the gravity component alone while dashed lines show the tree-level bispectrum 
including the primordial component for the local (upper right panel), equilateral (lower left panel) and orthogonal (lower right 
panel) models assuming the values of /nl corresponding to the 95% C.L. limits as determined by |Smith et al.^ i ,2009) and |Senatore| 
\et a/.|(2009l l from WMAP observations. The shaded area indicates the currently allowed region. 



models, for values of the respective parameters /nl corresponding to the current 95% C.L. limits ( Senatore et at. 
[2009; Smith et al. , 2009 ), with the shaded area indicating the allowed region. 

In addition, presents a specific dependence on triangle shapes, determined by gravitational instability and 



described by ( IV. 167 i at tree-level. The shape dependence of Bo, determined by the specific non-Gaussian model under 
consideration, is generically different. Such differences can be explicitly shown in plots of the reduced bispectrum, 
defined as 



B(kuk2,h) 



P(/ti)f(/t2) + 2perm. 



(IV. 170) 



which removes the redshift and scale dependencies of the gravity contribution, fig. 



16 



shows the reduced bispectrum 
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FIG. 16 Effect of the primordial component for different non-Gaussian models on the matter reduced bispectrum, as a function of 
the triangle shape. The continuous line shows the reduced bispectrum Q(k[ , fe, fcj) at tree-level in PT for Gaussian initial conditions 
at redshift z = 1 assuming k[ = 0.01 /iMpc"', = 1.5ki as a function of the angle 6 between ki and fc- Dashed lines show the 
reduced bispectrum including the primordial component for the local (upper left panel), equilateral (upper right panel), orthogonal 
(lower left panel) and folded (upper left panel) models. For the local, equilateral and orthogonal models we assume the values of 
/nl corresponding to the 95% C.L. limits as determined by |Smith et al. | ( [2009 ll and'Sena tore ef fl/.l ( |2009[ l from WMAP observations. 
The shaded area indicates the currently allowed region. For the folded model, for which no observational constraints are available, 
the values /nl = ±300 are considered. 



Q{k\,k2,k-i) at tree-level in perturbation theory, at z = 1 for - 0.01 /zMpc"', ^2 - l-5A:i as a function of the angle 
9 between ki and k2- In all panels, the continuous line represents the gravity-induced term which assumes larger 
values for nearly collapsed triangles, i.e. for ^ or n. This indicates that the probability of finding larger values 
for the matter density in triplets of points forming a squeezed or folded triangle is larger than for nearly equilateral 
triangles. This prediction is confirmed by the typical filamentary nature of the large-scale structure, evident from 
snapshots of N-body simulations or images of redshift surveys, since along these filaments it is easier to form collapsed 
triangles than equilateral ones. It should be stressed that the bispectrum is, in fact, the lowest order statistic sensitive 
to the three-dimensionality of structures and that these features are not captured by the information contained in 
the power spectrum alone. The effects of the primordial component on the matter bispectrum are shown by the 
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dashed lines which correspond, as in fig.[T5]to the 2-cr limits from CMB observations, in the case of the local {upper 
left panel), equilateral {upper right panel) and orthogonal {lower left panel) models while they correspond to the 
values /nl = ±300 in the folded case, for which no experimental bounds are available. Although the large scales 
k\ - 0.01 /iMpc ' and k2 - 0.015 /z Mpc ' and the relatively high redshift z - I have been chosen to enhance the 
effect of the non-Gaussian component, these triangles are not completely out of reach for future, large-volume surveys. 
Primordial non-Gaussianity modifies, in very specific ways, the shape dependence of the matter bispectmm produced 
by gravitational instability. 

While the dependence of the matter bispectrum on scale and redshift is responsible for the specific behavior of 
the skewness of the matter density field on the smoothing scale R and redshift, the sensitivity to the triangle shape is 
completely lost in analysis of the density higher-order moments. Instead, accurate measurements of the bispectrum, 
when achievable, offer in principle the possibility to disentangle the different contributions when triangles of different 
size and shape are included in the analysis. 

The matter bispectrum is not, unfortunately, a direct observable. While we will discuss later how the statistical 
properties of the matter distribution can be inferred from galaxy redshift surveys, we should mention that the shear 
field in weak lensing surveys is another observable directly related to the matter distribution. The observational 
consequences on the weak lensing bispectrum, of a primordial non-Gaussian component (of the local type) such as the 
one in (IV.168 1, have been explored in Takada & Jain (2004i. The authors find that, the primordial component alone 
{i.e. without contamination from the gravitational one) could be detected if > I50f^^^, assuming /,„n.v - 500 and a 

tomography over four redshift bins for a galaxy number density of tig = 100 arcmin"^. The large cosmic variance for 
low /"s makes difficult the detection of the primordial component, prominent instead at larger scales. As we will see 
in the next section, primordial non-Gaussianity has some effect on small scales as well, due to the nonUnear evolution 
of structures. 



2. Second-order corrections 

The simple prediction of ( |IV.168[ ) for the matter bispectrum is expected to be valid at the largest observable scales 
and at high-redshift, where nonlinear evolution is subdominant. Despite the fact that such conditions correspond as 
well to the regime where a detection of the initial component Bq is favored, the effects of non-Gaussian initial con- 
ditions can be significant even at smaller scales and at low redshift. Since these effects are the result of nonlinear 
gravitational evolution and non-Gaussian initial conditions, it is no longer possible to identify distinct contributions 



resulting from distinct sources of non-Gaussianity, as it is the case for the tree-level expression of (IV.168 i. Neverthe- 
less, it is possible to distinguish individual corrections in PT to the matter bispectrum depending exclusively on the 
initial power spectrum Pq, and therefore present as well for Gaussian initial conditions, and corrections depending in- 
stead on higher-order initial correlators, such as the initial bispectrum Bq and trispectrum Tq, which can be interpreted 
as small-scales effects due to non-Gaussian initial conditions. One-loop corrections in PT for Gaussian initial condi- 
tions have been studie d in | Scoccimarro| ( |1997l l, while the extension of these results to non-Gaussian initial conditions 



is studied in'Sefusatti|( |2009 



A comparison of these results with measurements of the matter bispectrum in N-body simulat ions ([Desjacques et al. 



20091 with non-Gaussian initial conditions of the local kind can be found in |Sefusatti et al. (20101. Fig.|17|shows 



the equilateral configurations of the matter bispectrum measured in N-body simulations together with predictions 
from perturbation theory at tree-level {dashed line) and one-loop {continuous line). In particular, the upper left panel 
considers B{k, k, k) for Gaussian initial conditions while the upper right panels shows the same quantity divided by 
the tree-level prediction in PT to highlight the small-scales non-linear behavior. The lower left and right panels show, 
respectively, the ratio and the difference between the matter bispectrum with an initial local component corresponding 
to /nl = 100 and the Gaussian case. The agreement between one-loop predictions and the simulations results is 
quite remarkable, while we notice that the tree-level prediction fails to accurately describe the effect of primordial 
non-Gaussianity already at relatively large scales. 

The significance of these relatively small corrections to individual configurations is to be considered in relation to 
the much larger number of configurations that can be measured as we include smaller and smaller scales and they could 
lead to a measurable effect when considered in terms of the cumulative signal to noise ratio. On the other hand, these 
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FIG. 17 Upper panels: equilateral configurations of the matter bispectrum measured in N-body simulations with Gaussian initial 
conditions (data points) and tree-level (dashed lines) and one-loop (continuous lines) predictions in perturbation theory. The right 
panel shows the ratio to the tree-level prediction with acoustic oscillations removed. Lower panels): ratio (left) and difference 
(right) between the matter bispectrum measured in realizations with local non-Gaussian initial conditions (/nl = 100) and the 
Gaussian case, compared with PT predictions. From ( |Sefusatti et a/.|[2010^ . 



effects loose in part the shape dependence of the original initial bispectrum and require an accurate model (perhaps 
beyond standard perturbation theory) and strong priors on the underlying cosmological parameters to be distinguished 
from the nonlinear, "Gaussian" component. A step in the direction of improved predictions is offered by the promising 
results of Renormalized Perturbation Theory (Bemardeau et at. |2008[ Crocce & Scoccimarro 2006a|b i and of the 
Renormalization Group approach ( [Matarrese & Pietroni 2007 [ Pietroni 2008 1. The extension of the latter to the 
case of non-Gaussian initial conditions has been recently considered in Bartolo et at. ( 2009| l, which studies specific 
predictions for the matter power spectrum and bispectrum. 



C. The Galaxy Bispectrum 



From the discussion above, we could expect that future, large-volume and high-redshift galaxy surveys will be able 
to directly detect a possible, large primordial component to the matter bispectrum by measurements of the galaxy 
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bispectmm, or at least provide constraints on the non-Gaussian parameters comparable to the constraints from mea- 
surements of the CMB bispectrum. Such an expectation is motivated by the simple observation that the number of 
Fourier modes available in a three-dimensional, ideal, all-sky galaxy survey is in principle much larger than the number 
of modes available in two-dimensional CMB maps. 

The galaxy distribution is, however, a less direct probe of the early Universe than the CMB temperature fluctuations. 
On top of the nonlinear evolution of structures and its contribution to higher-order correlation functions, one has to 
take into account the nonlinear nature of galaxy bias, itself responsible for additional non-Gaussianity. An analysis 
of the galaxy bispectrum should therefore be able to detect a small primordial component by separating it from these 
primary contributions. 

In this respect, an even more complex picture, due to additional and somehow unexpected effects of primordial 



non-Gaussianity on galaxy bias, has been emerging in the last couple of years, following the results of Dalai et al. 



( |2008| l. N-body simulations have shown, in fact, that nonlinear bias and an initial bispectrum are not two distinct 
sources of non-Gaussianity for the galaxy bispectrum, not even at large scales! Instead a local initial component 
can significantly affect the bias relation precisely at large scales, adding extra corrections. In the spirit of a review 
and since we do not have, at the time of writing, a satisfactory model of the galaxy bispectrum in presence of non- 
Gaussian initial conditions, in this section we will summarize earlier results, while in section [rV.C.5| we will present 
the recent developments that radically changed our understanding of the effects of local non-Gaussianity on the large- 
scale structure and finally comment, in section 1V.C.6 on some consequences for galaxy bispectrum measurements as 
far as current research provides. 



1 . The galaxy bispectrum and local bias 



Until recently, it was commonly assumed, even for non-Gaussian initial conditions, that the galaxy overdensity 
5g(x), defined in terms of the galaxy density ng{x) and its mean hg as 



(5„(x) = 



ng(x) - n„ 



(IV.171) 



can be expressed, at large-scales, as a local function of the matter density contrast, (J(x)'**, i.e. 

6g(x) = f[6(x)] . 



(IV. 172) 



Such a reasonable expectation is based on the fact that the physics of galaxy formation operates on much smaller 
scales, below the typical halo size, than those we are interested in. At large scales, where fluctuations are small, 
Sr < 1, we can consider the Taylor expansion, (Fry & Gaztanaga 1993 i 



1 



1 



6g{x) = bi d(x) +-b2 6\x) +- b3 6\x) + . 



(IV. 173) 



describing the bias relation between galaxy and matter in terms of a series of constant bias parameters, bj. This 
expansion allows for a consistent extension of the perturbative expressions for the matter correlators to the galaxy 
ones. In fact, from (IV.173 1 we can derive the galaxy three-point function in position space 



{6g{xi)6g(x2)6g{x3)) = b] {6(xi)6(x2)6(x3)) + b]b2 {6(xi)6(x2)5\x3)) + perm. + . . 
and the tree-level expression for the galaxy bispectrum given by 

Bgiki , k2, h) = b\ B"'Hki , k2, h) + b\b2 [PQ(ki)PQ(k2) + 2 perm.] , 



(IV. 174) 



(IV. 175) 



Properly speaking we should consider here the smoothed matter density contrast, that is 6r(x) = J c/^x' Wr(x - x')i5(x') with Wr a top-hat filter 
function. For simplicity, we implicitly assume a smooth density field, so that, for large enough filtering scale, e.g. R ~ 10 /i"' Mpc, matter 
perturbations are small, 6^1 
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where the second term on the r.h.s., proportional to the quadrat ic bias parameter b2, is of the same order of the gravity- 
induced contribution to the matter bispectrum B'^'^'', (IV. 167 1. Relying on this simple result, measurements of the 



galaxy bispectrum has been considered in the first place, m the context of Gaussian initial conditions, as a way to 
determine the bias parameters, and break the degeneracy between linear bias (^i) and the amphtude of matter fluctu- 
ations (e.g. cTg), otherwise affecting power spectrum measurements (Fry 1994!'Jeong & Komatsu 2009a Matarrese 
etar\\l 997; Nishim ichTeFar] [20071 [Scoccirnarro[ |2000a[ [Scoccimairo etal.\ \1999: Sef usatti et a/.||2006[|Sefusatti & 
Scoccimarro| |2005| l. In this respect, the corresponding reduced galaxy bispectrum is 



P^(^i)Pg(fe2) + 2 perm. bi b] 



(IV. 176) 



where Q is the reduced matter bispectrum (including a possible initial contribution) and the eff'ect of nonlinear bias 
is simply given by an additive constant term. As already mentioned, measurements of triangular configurations dif- 
ferent in shape and size allow to disentangle the different sources of non-Gaussianity and determine independently b\ 
and i>2, provided that accurate predictions for the matter bispectrum, from PT or N-body simulations, are available 



( Guo & Jing , 2009a b) and the effects of redshift distortions and the survey geometry are properly taken into account 
|Scoccimarron 2000a; S mith ef a/.||2008l ). 

In particular, if we allow the possibility of non-Gaussian initial conditions, then the matter bispectrum includes an 
initial contribution, so that we can rewrite (??) at tree-level explicitly as 



QT = jr [e/(/NL) + QT\ 



b_2 

b^ 



(IV. 177) 



and we can extend the analysis to obtain simultaneous constraints on the bias parameters and on the parameter deter- 
mining the amplitude of the primordial bispectrum, i.e. /nl- A first conservative estimate of the possibiUties offered 



by this method in measurements of the galaxy bispectrum in the 2dF Galaxy Redshift Survey (CoUess et al. 2001 



and in the Sloan Digital Sky Survey (SDSS York et al. 2000| l is given in | Verde et a/.| ( [2000) l as a simple extension of 
previous results for the bias alone ("Matarrese et al. 1997 1 suggesting that a primordial component could be detected 
for values of a local /nl of the order of lO^'-lO'*. As we will see in the next sections, a complete analysis of the galaxy 
bispectrum, including all measurable configurations can improve this estimate by more than an order of magnitude. 

Among the various observational issues in analyses of galaxy correlators, e.g. finite volume effects or completeness 
of the galaxy samples, we stress that particularly relevance has the problem of redshift distortions. Redshift distortions 



have in fact a significant impact on the shape dependence of the galaxy bispectrum, particularly at small scales ( Scoc 



cimarro 



2000a Scoccimarro et al. 1999 j l. A recent treatment of redshift distortions in bispectrum predictions (with 



Gaussian initial conditions) can be found in |Smith et a/.| ( |2008l l 



2. A bispectrum estimator 

In this section we define a simple estimator for the measurement of the galaxy bispectrum in N-body simulations as 
well as actual data. This allows us to derive an expression for the bispectrum variance and define a Fisher matrix for 
an analysis of the galaxy bispectrum in terms of the non-Gaussian (and bias) parameters. In the next section we will 
consider a proper likelihood analysis and the effects of the bispectrum covariance. Since what follows can be applied 
in general to bispectrum measurements we will consider, to simplify the notation, the case of the matter density field 
in Fourier space, described by the density contrast 6]^. We will point-out relevant differences in the application to the 
galaxy distribution. 

For a cubic box of volume V, a bispectrum estimator can be defined as ( Scoccimarro et al. 1998[ l 



B(h,k2,h) = {d^q2 fr<?3 5D(qi23) <Jq. \, , (IV.178) 

VB(ki,k2,k2) Jk, Jk, 

where Vf = k^^. - (Inf' jV is the volume of the fundamental cell and where each integration is defined over the bin 
qi e [ki - Ak/2,ki + Ak/2] centered at ki and of size Ak equal to a multiple of the fundamental frequency kf. The 
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Dirac delta function 6d(<Ii23) ensures that the wavenumbers qi, q2 and q3 indeed form a closed triangle, as imposed 
by translational invariance, while the normalization factor Vsiki , k2, kj,), given by 



VBiki,k2,kj) = Id^qi ld^q2 Id^qidDiqm) - kik2k3 Ak^ , 

J k] *Jk-> *Jk-t 



(IV. 179) 



represents the number of fundamental triangular configurations (given by the triplet qi, q2 and qs) that belong to the 
triangular configuration bin defined by the triangle sizes ki, k2 and kj, with uncertainty Ak. 

The lea ding contribution to th e bispectru m variance following from this estimator, in analogy with the power spec- 
trum case ( jFeldman et al. 1994 1, is given by Scoccimarro et al. ( 1998 1 

AB\ki,k2,h)^Vf- 



with the factor s 



T7 /7 7 1 \ '■tot {ki)P,o,(k2)P,o,(k3) , 

VB(ki,k2M) 

123 = 6, 2, 1 respectively for equilateral, isosceles and general triangles and where 

Ptotik) = P(k) + ^ ^ 



{2ny n 



(IV. 180) 



(IV.181) 



with the particle (or galaxy) number density n accounting for the shot noise contribution. In the case of a galaxy 
distribution, the matter power spectrum P{k) on the rh.s. should be replac e with th e galaxy p ower spectrum, expressed. 



at large-scales, by Pg{k) - b\P{k), under the local bias assumption of (|IV.173|. (|IV.180|l constitutes the Gaussian 



limit to the bispectrum variance, as it neglects higher-order corrections dependent on the three-, four- and six-point, 
connected, correlation functions. 



3. Fisher matrix forecasts 



In this section we consider simple forecasts for the constraints on the non-Gaussian parameters from measurements 
of the galaxy bispectrum in future redshift surveys. Specifically, we will consider a Fisher matrix for reduced galaxy 
bispectrum Qg in terms of the non-Gaussian parameter /nl and the linear and quadratic bias parameters bi and b2- 
These three parameters characterize the relative weight of the different non-Gaussian contributions to the galaxy bis- 
pectrum. Since the possibility to detect a primordial component relies on our ability to separate the three contributions, 
a robust result should, at least, involve a marginalization over bias. On the other hand, we will assume all cosmological 
parameters as known. This is in part justified by the weak dependence of the matter reduced bispectrum on cosmology 
discussed in the previous section. In this respect, it can be shown that the reduced bispectrum has the same signal to 
noise as the bispectrum. For a given triangular configurations, in fact. 



Qg{ki,k2,k-i) Bg{kx,k2,k'i) 



AQg{kuk2,ki) ABg{kuk2,ki) 



(IV. 182) 



since the variance of Q is dominated by the variance of B (see, for instance, Scoccimarro et al. 2004| l. 
The Fisher matrix can be written as 



triangles 



OQgdQg 1 
dpa dpp AQj 



(IV. 183) 



where the indeces a and /3 run over the parameters of interest /nl, bi and b2, while the reduced bispectrum variance, 
as mentioned above, can be expressed in first approximation as 



AQl(kuk2,k,) 



ABl{ki,k2,h) 



[Pg{h)Pg{k2) + 2 perm.]2 



(IV. 184) 



This expression, (see Scoccimarro et al. 2004 , corrects a typo in equation (A16) of Scoccimarro et al. 



1998 
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FIG. 1 8 Upper panels: predicted errors on galaxy bias parameters b\ (left) and Z?2 (right) as a function of the maximum wavenumber 
^max considered for the sum defining the Fisher matrix, l |IV. 185[ >. The analysis corresponds to an ideal geometry survey of volume 
V = lOh^^ Gpc' and a galaxy number density of Hg = 5 x 10"' Mpc"'. Dashed (red) lines assume a mean redshift of z = 1, 
while continuous (blue) lines assume z = 3. Both assume Gaussian initial conditions, i.e. /nl = 0. The vertical lines correspond to 
the value of k^.^^ determined as the inverse of the distance scale R defined by the condition (t(R,z) = 0.5. Lower panels: predicted 
errors on the non-Gaussian parameters f^- (left) and f^ (right), marginalized over the bias parameters, as a function of A:„,ax. From 
[Sefusatti & Komatsu|p007} . 



with AB^ given by (IV.I8O1. Notice that AQ^ depends on the linear bias parameter b^. The sum over the triangles 



configurations can be explicitly defined in terms of three sums over the wavenumber ki, ^2 and ^3 in steps of Ak, 



Z = Z E Z. 



triangles ki =kn^m ^'2=^m 



(IV. 185) 



with k*^.^^ - max(A-min, l^i - ^2!) to ensure that a close triangle can be formed and with ^^ax representing the minimal 
physical scale included in the analysis. Clearly, larger values of fc^j^x correspond to a much larger number of available 
configurations. For this reason, in fact, the cumulative signal to noise for the bispectrum, i.e. the sum of the signal 
to noise over all measurable configurations, grows more rapidly with fc,nax than it does for the power spectrum. On 
the other hand, we expect the primordial component to decrease significantly at small scales (high-fc). In practice, 
however, A;n,ax can be defined as the smallest scale at which we can trust our model for the galaxy bispectrum, in our 



case, the tree-level expression in (IV. 177 1. 

In fig. [18] the forecasted errors on bias parameters and non-Gaussian parameters as a function of k„ 



for an ideal 



geometry galaxy survey of volume V 
redshift z = 1 (dashed, red lines) and z 



10 h Gpc-' and a galaxy number density of = 5 x 10 ■'/i-'Mpc ^ at 
3 (continuous, blue lines). The negligible difi'erence between the results for 
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the non-Gaussian parameters at different redshift is a consequence of the fact that the signal to noise of the primordial 
component to the matter and galaxy bispectrum for a single triangular configuration, Bq/AB^, is, in our approximation, 
constant, both as a function of redshift and scale. This is not the case for the contributions due to gravitational 



instability and bias. It is clear that the choice of fe^ax significantly affects the final result. For instance, Sefusatti & 



Komatsu ( 2007 1 define k^^^, for a given survey, from as the inverse of the scale R given by the condition cr{R, z) - 0.5 



max 



to ensure that the tree-level predictions is applied within the mildly nonlinear range. Notice that the choice of k, 
depend on redshift, since at larger redshift we can expect a larger range of validity of perturbation theory predictions, 
both for matter and galaxy bispectrum. The dependence on the survey volume is simply given by - 1/VV. 

Scoccimarro et al. ( 2004| l, from a Fisher matrix analysis as the one described above, have shown that the 2dF and 

' ' ' ' f\oc ' ' ' ' " ' 



SDSS surveys should be able to probe values of /'°'^- < 100, assuming fcmax - 0.3/iMpc They suggested as well 



that an all-sky survey with a galaxy number density of ~ 3 x 10"^ Ir" Mpc"^ up to redshift z ~ 1 can probe values of 
f m.' of order unity. 

Sefusatti & Komatsu (2007) provided more specific predictions for a choice of planned and proposed high-redshift 



galaxy surveys, based on a similar Fisher approach, for the errors on non-Gaussian parameters both for the local 
and equilateral model. It is found that, for equilateral non-Gaussianity, the degeneracy between the non-Gaussian 
parameter /^^ and the bias parameters is severe, but it extends to unphysical regions of the bi-b2 plane and it can 
be severely reduced by introducing a correlation between linear and quadratic bias as the one predicted by the halo 
model. The marginalization over bias can be then replaced by a marginaUzation over the parameters of the Halo 
Occupation Distribution describing the galaxy population. Sefusatti & Komatsu ( 2QQ7\ finds that future large-volume 
surveys (V ~ 100 Gpc^ at z ~ 1,2), designed to accurately measure acoustic oscillations in the galaxy correlation 
function and thus map the late-time expansion of the Universe, should be able to probe ~ 4 and f^ ~ 20, 
i.e. values comparable to those expected from future CMB missions. At that time they constituted they best forecasts 
for constraints on /nl from large-scale structure measurements. These results implied, in particular, that if Planck will 
indeed detect primordial non-Gaussianity, a confirmation by large-scale structure observations will be required. 



4. Effects of covariance and current results 



The simple Fisher matrix analysis described in the previous section makes several approximations, starting with 
the assumption of an ideal geometry for the survey under consideration, and the Gaussian variance for the galaxy 
bispectrum configurations. In fact we can expect a proper treatment of the survey selection function and of the 
bispectrum covariance to have a significant impact on the estimation of the non-Gaussian (and bias) parameters. For 
instance, triangular bispectrum configurations at the largest scales probed by a realistic redshift survey (where the 
initial component should provide the largest corrections) are indeed highly correlated, because of the limited number 
of measurable Fourier modes. 



The issue of bispectrum covariance has been studied in ( Scoccimarro 2000a| 



et al. 2006[ Sefusatti & Scoccimarro 2005 | l. For instance, Scoccimarro et al. P004| l compare the Fisher matrix 



Scoccimarro et al. 2004 Sefusatti 



results for an ideal survey with a volume and galaxy number density similar to those of the main sample of the 
SDSS, with the predictions resulting from a likelihood analysis of the same survey, including the effects of survey 
geometry and covariance. Such analysis involves all measurable triangular configurations defined by wavenumbers ki, 
ki, ki < 0.3 h Mpc"\ with - 0.015 h Mpc"', resulting in a total number of triangle bins, Nj = 1015. The estimation 
of the corresponding, 1015 X 1015, bispectrum covariance matrix clearly represents a challenging computational 
problem as it cannot be determined from a relatively small number of N-body simulations. This work uses instead a 
code ( Scoccimarro 2000a| l implementing particle displacements as predicted by second order Lagrangian perturbation 



theory (2LPT, see, for instance, Bernardeau et al. 2002 and references therein) to produce 6,000 realizations of the 
density field. Such large number of realizations is in fact necessary for an accurate determination of the covariance 
matrix. In addition, the 2LPT results, including particle velocities, allow for exact redshift distortions. Each mock 
catalog, in redshift space, is then weighted according to the FKP procedure ( [Feldman et al.\ |1994[ |Matarrese et al.\ 
1997 Scoccimarro 2000a|l to take into ac count the SDSS selection function. The same covariance matrix is compared 



to analytic expressions in Sefusatti et al. ( 2006[ l 



Given a proper estimate of the covariance matrix, a likelihood function for the reduced bispectrum 2„ can be defined 
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FIG. 19 95% confidence limits on /nl from the PSCz galaxy bispectrum after marginalization over bias parameters, as a function 
of the number of eigenmodes included in the likelihood analysis. From |Scoccimarro et a/.| ( (2004[ ). 



in terms of the normalized bispectrum eigenmodes q„ that diagonalize it ( Scoccimarro 2000a i. These can be expressed 

as 



— ^ ymn 

m=l 



Qm-Q 



(IV. 186) 



where = (Qm), ^Ql, = {(Qm - Qmf) and their signal to noise is given by 



\N)n 



1 



Nt 

^ \ ymn 



Qm 



(IV. 187) 



where A„ represents the eigenvalue for q„, with {q„ q,„) - 6„„,. The eigenmodes presenting the largest signal to noise 
can be easily interpreted by considering how they weight different bispectrum configurations. In fact, the largest signal 
to noise corresponds to an eigenmode defined by a nearly equal weighting of all triangles, and it therefore represents the 
overall bispectrum amplitude. The next eigenmode weights instead with opposite sign triangles close to the equilateral 
shape and nearly coUinear triangle. Each eigenmode represents in fact a fraction of the information contained in the 
bispectrum confi g uration s, and a crucial role in th is respect is played by the shape and scale dependence. To illustrate 
this point, fig. 



19 



(from Scoccimarro et a/. , 2004) shows the 95% C.L. limits on f^^- from the likelihood analysis of 
the IRAS PSCz catalog (Saunders et al. , 2000) as a function of the number of eigenmodes included. 

Although the diagonalization of the covariance matrix does not ensure the exact independence of the eigenmodes, 
which can still present non- vanishing higher-order correlations, it has been shown that this is nevertheless a reasonable 



assumption in practice ( Scoccimarro 2000a I . This allows us to write a likelihood function for the non-Gaussian and 
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bias parameters, denoted generically as pa, in terms of the product of the probability distribution functions P„(x) for 
each individual eigenmode, that is 



mPa]) 'X n PnlqniiPa})]. 



(IV. 188) 



n=l 



The probability distributions P„(x), which can be determined from the mock catalogs, are not expected in general 
to be Gaussian, although this can be in fact a good first-order approximation in the case of the SDSS main sample 
( |Scoccimarro et a/."||2004| l. 

A direct implementation of this kind of analysis, taking into account all measurable bispectrum configurations and 
their covariance, has been performed by Scoccimarro et al. ( 200 1 [) for different IRAS catalogs ( ,Efst athiou et al. 1990 



Fisher et aL\ |1995[ [Strauss et aL\ |1992| l and by |Feldman et aZ.| ( |200T) for the IRAS PSCz catal og (Sau nders et al. 



2000T considering the case of the x'^ model of primordial non-Gaussianity (|Scoccimarro| 



]2004) ) 



derives the limit 1/,^°'; | 



2000b I. Scoccimarro et al. 



< 1800 at 95% C.L. for the bispectrum measured in the PSCz catalog. 



Along this lines. 



Scoccimarro et al. 



be obtained from measurements of t 



(2004 1 also studied the constraints on f}?^- for local non-Gaussianity, that could 



le galaxy bispectrum in the SDSS main sample, including the effects of the 



survey geometry and bispectrum covariance, forecasting the 1-cr error ^f^^' - 150, after marginalization over the 
bias parameters. This work compared this more realistic estimate of the predicted errors on from the likelihood 
analysis of the SDSS bispectrum to the Fisher matrix forecast for an ideal geometry of nearly the same volume and 
galaxy density finding a worsening of a factor of 4-5. They point-out, however, that the realistic errors, which are an 
estimate from the north part of SDSS alone, should be taken as a an upper bound to the results actually achievable 
because of the FKP weighting scheme, not optimal at the largest scales where the primordial component is the largest 
and because of the fact that extra signal can be found as well in open configurations, not considered there, due to the 



broken translation invariance. We might add, based on the results of section IV.B.2 that nonlinear corrections present 
for non-Gaussian initial conditions might increase the overall signal due to a non-zero /nl, 
where a large number of triangular configurations can be measured. 



, particularly on small scales 

At this point we should remind the reader that all the results discussed so far on the galaxy bispectrum and its 



significance for constraining primordial non-Gaussianity, assume the expression (IV.177i to be a reliable prediction. 
As we shall se in the remainder of this section, this is not the case, as additional effects of non-Gaussian initial 
conditions have to be taken into account. Nevertheless, the primordial component, whose direct detection has been the 
main target of the earlier works discussed above, is still expected to provide a contribution to the galaxy bispectrum, 
and there are good reasons to believe that these results can be still interpreted as a "conservative estimate" of the 
possibilities offered by bispectrum measurements in the large-scale structure to test the Gaussianity of the initial 
conditions. 



5. Primordial non-Gaussianity and non-local Galaxy Bias 



The constraints and forecasts discussed so far in this section are based on the tree-level expression for the galaxy 
bispectrum, (IV.177i, derived under the assumption of local bias, (IV.173i. As anticipated, our understanding of 



galaxy bias in presence of primordial non-Gaussianity radically changed in the last two years, after Dalai et al. ( 2008 1 
presented measurements of the halo power spectrum in simulations with non-Gaussian initial conditions of the local 



kind showing the presence of large corrections at large scales, not captured by the local bias prescription! Fig. 20 
shows the matter-halo cross-power spectrum for different values of /nl from these simulations, where the unexpected 
effect of non-Gaussianity at large scales is evident. 



Local bias, (IV. 173 i), in fact, implies a leading contribution to the galaxy (or halo) power spectrum of the simple 
form 



PAk)^b\PQ{k), 



(IV 189) 



with no dependence on /nl, while the simulations results of Dalai et al. (12008'), later confirmed by Desjacques et al. 
( 2009] !; Grossi et al. ( 2009[ l; Pillepich et al. (|2008j), are consistent with a scale-dependent correction to the linear bias 
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FIG. 20 Matter-halo cross power spectrum measured in simulations with local non-Gaussian initial conditions for different values 
of /nl. From |Dalal e7aIH2b"08^ . 



of the form 



Pg(k) = [bi+Abi(k)fPoik), 



(IV. 190) 



with 



Abiik) = 3/nl(^i - m 



O 

k^T(k)D(z) ' 



(IV.191) 



where 6c - 1.68 is the linear critical density for spherical collapse, extrapolated at z = 0. Such coiTection therefore 
increases with the scale, with redshift via the growth factor D(z) and with the non-Gaussian parameter /nl, and 
vanishes for unbiased populations (bi - I). 
A theoretical interpretation, based on the peak-background split (Cole & Kaiser 1989|l has been assumed by Af- 



shordi & Tolley| ( [2008l l; 
different derivation, by 



Dalai et al 



( 2008|l ; Giannantonio & Porciani ( 2009 1; Slosarefa/. ( 2008] !, and, with a somehow 
(2008||. According to these works, the local relation between the galaxy density 



McDonald] ^ 

and the matter density in ( IV.173| l is modified, in presence of local primordial non-Gaussianity, to include an explicit 



dependence on the primordial curvature perturbation, <i>, i.e. 

1 



5„(x) = b,5{x) + ci(/nl)1>(x) + -b26\x) + C2(/nl)5(x)<D(x) + , 



(IV. 192) 



The galaxy two-point function, will be given by the following perturbative expansion, 

{6g(Xi)6g(x2)) = bl{6(xi)6(x2)) + bici(fy0{6(xi)<i>(x2)) + perm. 



(IV.193) 



where the second term can be rewritten as the scale-dependent of the linear bias pa rameter o f Eq. (|IV. 191|l, with the 
l/k^ behavior resulting from the relation between and Ok given by M{k) ~ k^. Giannantonio & Porciani, (2009( ) 



describes, in fact, the galaxy distribution as multivariate distribution, although the matter density 6 and the curvature 
O are not two independent random field, but they are related by Poisson equation, (IV.161 1. It should be noted, that 



no derivation of a similar effect due to a different kind of primordial non-Gaussianity (if feasible) has been, so far. 
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proposed. The derivations presented in the works cited above, in fact, all rely on the relatively simple expression 



defining the local model (II.6 1 while their generalization to a model defined by a generic initial bispectrum is a quite 
challenging problem. 

Following the results of Dalai et al. (2008 ), moreover, an apparently different explanation, resulting in fact in a very 
similar but distinct effect on the galaxy power spectrum has been proposed by Matarrese & Verde ( 2008| l and Taruya 



et al. (20081. Taruya et al. ( 2008| l, starting from the local bias prescription of (IV.173i, point-out that the next-to- 
leading order correction to the galaxy two-point function in presence of local primordial non-Gaussianity, represents, 



in fact, a large correction, identical up to a constant factor, in the large-scale limit, to the bias correction of (IV. 191 
The perturbative expression for the galaxy power spectrum is given by 



Pgik) ^ b\P{k) + bxbi d^q B(k, q, |k - q|). 



(IV. 194) 



where the second term, proportional to the quadratic bias parameter b2 and dependent on the matter bispectrum B, 
corresponds to the lowest order, one-loop correction. Remarkably, for local non-Gaussianity, in the limit k ^ Q, 
such correction presents the same scale and redshift dependence, and, for massive halos or highly biased populations 
{bi » 1), even the same amplitude, as the one resulting from (IV.191 1. The expression, however, can be applied 



to any model of primordial non-Gaussianity, given the appropriate initial matter bispectrum (see, for instance, | Verde| 
& Matarrese 2009| l. In the case of equilateral non-Gaussianity, the correction is almost negligible, while local non- 
Gaussianity appears t o be a limiting case le ading to a particularly significant effect. The same correction has been 
considered already by Scoccimarro ( 2000b i in the context of initial conditions, where it leads to a redefinition of 



the bias parameters, with no additional scale-dependence. 

Matarrese & Verde (2008i presents a different derivation of an expression similar to the one of ( |IV.194 1, based on 
earlier works on the density peak correlation function (Grinstein & Wise 1986^ Matarrese et al. |1986) . In this case, 
a specific prediction for the bias parameters, valid however only in the high density threshold Umit, is included. It 
is interesting to notice that the possibility of large-scale effects on the correlations of biased distributions has been 
explicitly pointed-out by Grinstein & Wise ([1986 ), although without further study. 

The two distinct corrections to the galaxy power spectrum, one corresponding to the modified bias relation of 



(IV. 192 1, the other to the perturbative correction due to nonlinear bias of (IV.194i, have been studied in a compre- 
hensive framework recently by |Giannantonio & Porciani| |2009 ), where the authors suggest that the effect measured 
in N-body simulations is mainly due to the multivariate nature of the galaxy distribution with local primordial non- 



Gaussianity, rather than the effect of nonlinear bias (IV.194i. In addition, Desjacques et al. (20091 pointed-out that 
even the galaxy bias parameters bi, related in the framework of the halo model to the halo bias parameters bh,i{M) 
for halo populations of mass M, present a dependence of /nl due to the effects of non-Gaussianity on the halo mass 
function. The picture that has been emerging in the last years is therefore quite complex and it should be stressed that 
a wide consensus in the community on a well defined model, even for the galaxy power spectrum, is still lacking. For 
instance, a discrepancy of the order of a 10% between predictions and simulations results, did not find yet a unique in- 
terpretation (see discussions in Deyacques'eFar 2009 Giannantonio & Porciani 2009 ; Grossi et al. 2009 |Maggiore] 



& Riotto 2009 , Pillepich et al. 



20081. 



This rather surprising effect of local non-Gaussianity on the bias relation leads, remarkably, to the possibility of 
placing limits on f}°^- from current large-scale structure observations, already comparable to limits from the CMB! 



( Afshordi & Tolley| |2008[ (Slosar et al. 2008j). Specifically, Slosar et al. (2008) derived from measurements of the 
cross-correlation of several large-scale structure data-sets and the CMB i^o et al.^(10Q%) the 2-cr constraints 



-29</^-<70, 



(IV. 195) 



leading to an marginal improvement of the WMAP results. Encouraging predictions for the constraints that can be 
derived in future spectroscopic as well as photometric redshift surveys can be found in Carbone et al. ( 2008[l. A f air 
comparison between these forecasts and those derived for the galaxy bispectrum in Sefusatti & Komatsu ( 2007[ l is 
clearly not possible as the latter do not include the effect on the bias relation discussed above. Two observations, 
however, are in order. In the first place, these effects on the galaxy power spectrum are specific of the local model 
of non-Gaussianity, while the galaxy bispectrum is in principle sensitive to any initial component Bq- In the second 
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FIG. 21 Large-scale contributions to the galaxy bispectrum due to primordial non-Gaussianity of the local (left panel) and equilat- 
eral (right panel) type described as one-loop corrections assuming a local bias prescription. Thin lines correspond to the contribu- 
tions for Gaussian initial conditions. From |Sefusatti|p009| l, see the reference for further details. 



place, robust results can be obtained from galaxy power spectrum measurements at large scales in photometric sur- 
veys. The degradation of the information that can be extracted from bispectrum measurements in photometric surveys 
with respect to spectroscopic ones is still to be properly studied. The impact of photometric errors on the accurate 
determination of the bispectrum dependence on the triangle shape can in fact be significant. 



6. The Galaxy Bispectrum after Dalai et al. (2008) 



First steps in the direction of an extension of the results discussed above to the galaxy bispectrum have been taken in 
Jeong & Komatsu (2009b I and Sefusatti (2009 1. Specifically, Jeong & Komatsu ( 2009b| l considered an expression for 



the high-peak three-point function derived in Matarrese et al. 



] 1986| l, analogous to the one for the two-point function 



studied by Matarrese & Verde (2008 ), and applied it to the case of local non-Gaussianity. Sefusatti (2009 ) considered 



instead the perturbative approach of Taruya et a/., ( ,2008, ) based on the local bias expansion of ( |IV.173| l, and applied it 
to local and equilateral non-Gaussianity. 

These works show that the galaxy bispectrum is expected to be sensitive to both the initial matter bispectrum Bq, as 
well as to the initial matter trispectrum Tq, by means of a contribution analogous to ( |IV.194[ ) and given by 



Br, ^ b\B(kuk2,h) + 



b\b2 



d^q T(ki,k2,q, |k3 -q|). 



(IV. 196) 



which represents a large correction at large scales, with an asymptotic behavior characterized by an extra 1 /A: factor 



Sefusatti 



(20091 



with respect to the primordial matter bispectrum component, Bq and a dependence on f^^. In addition, 
points-out that, unlike the power spectrum, large-scale corrections due to nonlinear bias are present as well for equi 
lateral non-Gaussianity (and virtually for any non-pathological form of the primordial bispectrum and trispectrum). 
Fig.|2T]shows the one-loop corrections to the galaxy bispectrum due to nonlinear bias and primordial non-Gaussianity 



under the assumption of local bias ("Sefusatti 2009 1. The left panel assumes local non-Gaussianity including a non- 



zero initial bispectrum and trispectrum, while the right panel assumes a non-zero initial bispectrum of the equilateral 
type. Thin lines coiTespond to Gaussian initial conditions. The black continuous line represents the matter bispectrum 
and therefore the first term on the r.h.s. of ( |IV. 196| l, while the blue dashed line coiTespond to the second term. Notice 
that at next-to-leading order in PT, the matter bispectrum T depends on the initial trispectrum as well as the initial 
bispectrum Bq, so that an effect is present also for equilateral non-Gaussianity where the figure assumes Tq = 0. 
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FIG. 22 Measurements of a set of triangular configurations of tire halo bispectrum in N-body simulations with local non-Gaussian 
initial conditions, as a function of the non-Gaussian parameter /nl- Large values of the parameter a correspond to m ore squeezed 
configurations. In the upper panels, the dependence of the halo bispectrum on f^^ is evident. From Nishimichi et al. ^2009^ . 



It should be noted, however, that these results ignore, at least for local non-Gaussianity, the modified bias relation 
of (IV. 192 1 (see, in this respect, some comments in [Giannantonio & Porciani 2009, ), and do not provide reliable 
predictions for the constant bias parameters. Furthermore, they have not been properly tested against measurements 
of the halo bispectrum in numerical simulations. The only work, at the time of writing, in this direction (Nishimichi] 
et al. 2009 ) shows, however that the dependence of the halo bispectrum on /nl is roughly consistent with the functional 



form resulting from the prediction of ( IV. 196 1. The authors attempt as well, using a simple fit to their measurements, a 
preliminary forecast analysis for a future large-volume (100 Z;"^ Gpc^), high-redshift survey, finding a detectable value 
of /nl - 20, using a very limited number of configurations. Fig.|22|from|Nishimichi et al.U2009\ shows measurements 



72 



of a set of triangular configurations of the halo bispectrum in simulations with local non-Gaussian initial conditions, 
as a function of the non-Gaussian parameter /nl, where the dependence of the halo bispectrum on f^-^ is evident. 

A simple but reasonable expectation would be that the inclusion of the effects of primordial non-Gaussianity on 
galaxy bias will improve the results of Sefusatti & Komatsu (2007 ), which are based on the detectability of the 
primordial component alone. Our understanding of these phenomena is, however, evolving rapidly in these days, 
and these notes on recent developments are likely to become outdated relatively soon. 



D. Running non-Gaussianity 

1 . The case of a scale-dependent /nl 

DBI models of inflation predict, as we have seen, a primordial curvature bispectrum very close to the equilateral 
model in its shape dependence. An additional but quite ge neric feature of these models is given by a significant 



departure from the hierarchical scaling B^{ k,k,k) ~ P%(k ) (Ch 



en '2005 



Chen et al. 



2007 



Khoury & Piazzaj |2008[ [Renaux-Petel , 2009,|Shandera & Tye„2006.). More recently, this po ssi 



Chen & Wang 



2009 



5ility has been explorec 



aswellinmodelsoflocalnon-Gaussianity (Byrnes ef a/. 2008 2009a|b Byrnes & Tasinato 2009[ Kumar et al. 2009 1. 

Under a phenomenological point of view, this extra scale-dependence can be described by a running /nl(^), or, 
more properly, in terms of an amplitude parameter /nl and a running parameter mnGj defined by 



(j^ \ hno 



(IV. 197) 



where kp is a properly chosen pivot scale, while K{ki,k2, ^3) = (ki + k2 + kj)/3 defines an overall scale characteristic 
of the triangular configuration on which B(^(k\,k2, kj) depends. In other terms, the /nl(^) defined above replaces the 
constant /nl in the definitions of the local and equilateral bispectra effectively introducing an extra dependence on 
scale. 

Observational consequences of a running /nl(^) have been explored in Lo Verde et a/. | ([2008) and Sef usatti et al. 
(j2009j), while in Taruya et al. ( 2008j ) this effect is included in the prediction for one-loop corrections to the matter and 
galaxy power spectrum. 

|Lo Verde et a l. (lOOW) provided an analysis of the possibility of constraining the running parameter hng by com- 
bining current limits from the CMB on the amplitude parameter /nl at the pivot scale kp = 0.04 Mpc"' with future 
measurements of cluster abundance. The effect of a «ng significantly different from 1, can result in a much larger (or 
smaller) amount of non-Gaussianity on the smaller scales relevant for the cluster mass function. Fig. 23 (left panel, 



from |Lo Verde et a/. 2008) illustrates the difference in the range of scales probed by different observables. Focusing 
in particular on the equilateral model for the curvature bispectrum, this work assumes the amplitude of /nl(^) to be 
constrained by the CMB bispectrum at the pivot point scale kp and derives the expected constraints on its running by 
considering the effective amplitude of /nl(^) at the smaller scales (k ~ 0.3 - 0.6/zMpc"') probed by cluster surveys. 
For an all-sky cluster survey up to redshift Zmax =13 they find the 1-cr constraints, marginalized over Q.„, erg and h, 
assuming the fiducial values /nl - 38 and hng - 0, Ahng - 2 with a Planck prior A/nl(^ - kp) = 40. Their analysis, 
however, does not include the simultaneous limits that measurement of the CMB bispectrum alone is expected to 
provide on both the amplitude /nl and running «ng. 



2. Running non-Gaussianity and bispectrum measurements 



Sefusatti et al. ( 2009| ) performs a Fisher matrix analysis of the CMB bispectrum to obtain the sensitivity of this 



observable to the running of /nl(^)- The results in the case of local non-Gaussianity, assuming the same pivot kp = 0.04 
and marginalizing over the amplitude /n°l% are the 1-cr uncertainties of A«ng ^ 0.68 (50//^°^" ) fsk^ WMAP and 
Ahng - 0.1 (50//^^ ) f^i;\!^ for Planck, where f^^- stands for the fiducial value of the amplitude parameter. In the 
case of equilateral non-Gaussianity, we have Ahng - 1-1 (100//^^)/''^^ for WMAP and Ahng - 0.3 (lOO//^?;) / ,'/^ 
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FIG. 23 Left panel: Range of scales probed by different observables compared with /nl(^) for different values of the running 
parameter h^q. From |Lo Verde et a/.| ( |2008) . Right panel: predictions of DBI models showing the peculiar relation between the 
amplitude /^^ and the running «ng for different values of the parameters of the inflaton Lagrangian. The figure incl udes the Fisher 
matrix forecasts for combined CMB and galaxy bispectrum measurements assuming a fiducial = -50. From 
j2009[ l, see the reference for further details. 



Sefusatti et al. 



for Planck. Since it is always possible, given the observable of interest (e.g. the CMB bispectrum for a specific 
experiment) and the non-Gaussian model, to choose the pivot point in such a way to remove any degeneracy between 
the amplitude and the running parameters, a measurement of the running parameter comes at no cost with respect to 
the determination of the /nl(^ = kp). Notice, however, that, for reasons related to the numerical implementation of the 
CMB estimator, |Sefusatti et al. ( 2009| l assumes for the overall scale representative of a given triangular configuration, 
the geometric mean of the three wavenumber, i.e. K = (^1^:2^3)'^^. While the difference with the more physically 
motivated definition in terms of the arithmetic mean K = (ki + k2 + ki,)/3 is very small for equilateral non-Gaussianity, 
in the local model this is not the case. 

[Sefusatti et aL\ ( |2009| l considers as well the Fisher matrix from large-scale structure information, and specifically 
the galaxy power spectrum (including the eff'ect on halo bias) for the local model and the galaxy bispectrum, but in 
terms of the simple description of (IV. 177 1, therefore excluding halo bias effects. This different choice of observables 
with respect to the model of primordial non-Gaussianity assumes a negligible effect of equilateral non-Gaussianity on 
the galaxy power spectrum (still to be confirmed by N-body simulations). It is shown, in particular, that future galaxy 
redshift surveys can significantly improve CMB results. Fig. 23 (right panel) shows the contours plots for the 1-cr 
uncertainties resulting for a joint Fisher matrix analysis of CMB and large-scale structure information. The expected 
limits are plotted against the predictions for the relation between the amplitude f^^- and the running «ng from DBI 
inflationary models. It is interesting to notice how these models predict a stronger running for smaller values of the 
amplitude parameter. In this respect, constraining the value of n^c can place additional limits on the parameters of the 
inflaton Lagrangian. 
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V. CONCLUSIONS 

Weakly non-Gaussian initial conditions are defined, in most of the relevant inflationary models, by a non-vanishing 
bispectrum for the primordial curvature perturbations. The most direct observables of this primordial density correlator 
are, naturally, the bispectrum of the temperature fluctuations in the CMB and the bispectrum of the mass distribution 
at large scales as probed by galaxy surveys. In this review we presented an overview of the problems, results and 
expectations connected with the detection of (or constraints on) primordial non-Gaussianity specifically in bispectrum 
measurements of the CMB and LSS. 

The CMB is an ideal observable for tests of primordial NG because temperature and polarization anisotropies can be 
described in the linear regime of cosmological perturbations. The statistical properties of the primordial curvature field 
are thus directly reflected in the pattern of CMB fluctuations. As we have seen, tests of primordial NG are formulated 
in terms of the estimation of the bispectrum amplitude /nl for each of the shapes predicted by different inflationary 
models. It was originally shown in the literature that a maximum likelihood estimator of the bispectrum optimally 
extracts all the /nl information from a CMB map. Extracting the primordial non-linear parameter from the bispectrum 
has subsequently become the standard way to test primordial NG in the CMB. The best /nl measurements to date come 
from analysis of the WMAP datasets, and roughly constrain the primordial bispectrum amplitude to be < 100 for the 
local, equilateral and orthogonal shapes. Despite already being very stringent (the NG part of the CMB temperature 
anisotropies is constrained at the level of 10"^ of the total fluctuation), these bounds are still far from the typical order 
of magnitude of primordial NG predicted by most inflationary models. As we have seen. Fisher matrix forecasts 
show that future results from the Planck satellite (whose release date is predicted to be in 2012) will improve previous 
WMAP constraints by roughly one order of magnitude, thus impacting the range of some theoretical predictions. 
This significant improvement is due to the better sensitivity of Planck to many more bispectrum configurations in the 
analysis, and the possibly of exploiting both temperature and polarization datasets. Another important limitation on 
current constraints is that inflationary predictions encompass more shapes than those that have been constrained so 
far. The reason why many shapes remain to be constrained is that they cannot be written as a separable product of 
one-dimensional functions of a single wavenumber Separability, as we have seen, is a crucial property since it makes 
the actual analysis computationally affordable in terms of CPU time. We have reviewed recent work showing that 
this limitation can also be overcome in future analysis by means of a fully general, and mathematically well defined, 
eigenmode expansion of the bispectrum shape. Thanks to this, and in light of the significant improvement in sensitivity 
provided by Planck, better and more general CMB constraints on primordial NG models will be available in the near 
future. One caveat is that the high precision of the forthcoming CMB datasets makes them much more sensitive to 
other spurious (i.e. non primordial) sources of NG, that could bias the /nl estimate. Achieving an accurate control on 
these contaminants is clearly a crucial goal for future analysis. As we have seen, much work is being done in order 
to predict, detect and isolate non-primordial NG eff'ects, but some issues still have to be addressed. In particular a 
complete prediction of the total bispectrum generated by second order cosmological perturbations is not yet available, 
although a number of effects have been studied in detail. Accurate characterization of NG from diffuse foreground 
residuals is another important issue that will require further investigation. 

For large-scale structure, many aspects of the general CMB scenario outlined above change, as should be evident 



from a comparison of the discussions in sections III and IV In the first place, we cannot rely on a direct relation 



between the observed galaxy bispectrum and the primordial curvature bispectrum predicted by inflationary models. As 
we have seen, a small departure from Gaussian initial conditions should result in a correction to the galaxy bispectrum 
induced by gravitational instability and nonlinear bias, constituting the dominant contributions. The nature of this 
correction is a complex problem in its own right, since it is due to the linearly evolved initial matter bispectrum as well 
as to the effects of primordial non-Gaussianity on the galaxy bias relation. Such effects are still under investigations 
and we do not have, to date, an accurate theoretical model. On the other hand, early results from galaxy power 
spectrum measurements are very encouraging, albeit restricted at present to the local non-Gaussian model. Current 
data sets already appear to be able to confirm and improve CMB results. In this respect, it is evident that the ultimate 
goal is the implementation of a complete large-scale structure analysis in terms of all measurable correlators, including 
power spectrum, bispectrum and beyond, i.e. an analysis that fully reflects the non-Gaussian nature of the mass and 
galaxy distributions even on large scales. 

There are several issues which remain to be resolved, for which we can identify three main categories. First, we need 
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to develop a robust model for the galaxy correlators accurately accounting for small-scale nonlinearities for both the 
matter and galaxy density fields, as well as in the presence of non-Gaussian initial conditions; this also must account 
describe non-localities in the bias relation. In this review we have briefly summarized the state of the art, noting that 
our understanding of these phenomena is evolving rapidly. Secondly, once a reliable model is available, it will be 
necessary to develop the machinery that will allow us, in the event of a future detection, to properly identify the effects 
of diff'erent models and their bispectrum shapes. In this respect, the CMB results presented in section III provide an 
important benchmark. Finally, observational problems connected with redshift surveys such as the effects of redshift 
distortions and/or photometric errors, survey selection function, completeness, etc., will have to be addressed. We have 
not discussed these issues here as they are generic to all large-scale structure experiments, but they clearly represent 
a major challenge for the exploitation of future data-sets. Both the first and the last point are crucial for virtually all 
the science goals of future ground-based or satellite surveys, particularly dark energy studies. Although only partial 
results have been obtained so far, there is every indication that characterising non-Gaussianity in future galaxy surveys 
will result in a significant test of the initial conditions of the Universe. 

To summarize, sufficient experimental sensitivity has been reached recently in CMB experiments (namely WMAP) 
to allow for meaningful constraints on the non-linear parameter /nl for several different families of models. These 
results are already arguably the most stringent quantitative test of the predictions of standard inflation. However, much 
tighter constraints on a broader range of models are expected from the future Planck data release. Thus a dramatic 
confrontation is set to continue between the de facto standard model of inflation and observational datasets from 
both the CMB and large-scale structure. Tests of primordial non-Gaussianity are rapidly becoming one of the most 
effective and promising approaches for gleaning important information about the physical processes that generated the 
primordial cosmological perturbations. 
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Appendix A: Basics of estimation theory 

If a random variable x is characterized by a Probability Distribution Function (PDF) /?(x|/l) dependent on a parameter 
A then an estimator for /I is a function £(x) used to infer the value of the parameter. If a given data set {x"*') is drawn 
from the distribution p{ii. A), then A - £(x°''') is the estimate of the parameter A from the given observations. Since 
fi is a function of a random variable, it is itself a random variable. In the literature a random variable obtained as a 
function of another set of random variables is often referred to as a statistic. 

A general property usually required when building an estimator is its unbiasedness. An estimator for a parameter A 
is unbiased if its average value is equal to the true value of the parameter: 

(A) = A. (A.l) 

The standard deviation is generally used to determine the error bars on A i.e. 




(A.2) 



where (.) denotes statistical average and cr^ is the variance of the inferred parameter When we measure a parameter 
A from a set of observations drawn from the PDF p{ii\A), we clearly would like our estimate not only to be unbiased, 
but also to have as small error bars as possible. In other words, among all the possible unbiased estimators of A that 
can be built, we look for the one that minimizes cr i defined in \A.2) . If such an estimator exists, it is called an optimal 
estimator. 



76 



In this context a crucial role is played by the Fisher information matrix, defined as 

The Fisher matrix appears in an important theorem, known as the Cramer-Rao inequality, stating that /or any unbiased 
estimator of A 

o-A > -^=- (A.4) 



This theorem is then placing a lower bound on the error bars that can be attained when estimating a given parameter 
from a given set of observations. No matter which estimator is used, the smallest attainable error bars will be given by 



the square root of the inverse of the Fisher matrix. For a demonstration of this crucial result see e.g. Kendall & Stuart 



1979|) or, in relation to the CMB bispectrum, Babich (2005 ). It is then clear that the best estimator of a parameter is 



an unbiased estimator saturating the Rao-Cramer bound. If such an estimator is found, then it is impossible to obtain 
a better estimate using any other statistic. The question then becomes if, for a given the PDF p(x\A), an estimator 
saturating the Rao-Cramer bound exists. 

It can be shown that a necessary and sufficient condition for an estimator £(x) of a parameter A to be optimal is the 
following: 

'-^-FAAm-A), (A.5) 
oA 

where F is the Fisher information matrix just introduced above. 

Another crucial quantity in estimation theory is the so called maximum-likelihood estimator. In a maximum- 
likelihood (ML) approach we take the observed data set x"*' as fixed and we estimate A as the parameter that maximize 
the probability (likelihood) to observe the given data. In formulae, the ML-estimate of A is the value A that satisfies: 

In this context the PDF p{x\A) is often denoted as the likelihood function and indicated as X.{x, A). Two powerful 
theorems involving the likelihood have been proven: 

1. If there is an optimal unbiased estimator (i.e. an unbiased estimator saturating the Rao-Cramer bound) then it is 
the maximum-likelihood estimator or a function of it. 



2. The maximum likelihood estimator is asymptotically optimal, i.e. it saturates the Rao-Cramer bound when 



— > oo, N being the number of repeated observations in our data set x"^,\ . . . , x"** 



These two theorems answer our initial question about the best estimator choice. The first theorem basically states that 
if a best method exist, then the ML-estimator is that method. Note that this result follows naturally from the optimality 



condition (A.5 1 introduced above. The second theorem says that for very large data sets the ML-estimator is the best 
method, i.e. the one saturating the Rao-Cramer bound. In other words, when dealing with the practical problem of 
estimating a parameter from a given data set, we should in theory always choose a ML-likelihood approach. However 
in practice this is not always possible: for example, the PDF p{x\A) might be too difficult to calculate or sample 



numerically, or the ML condition ( A. 6 1 (generally a complicated non-linear equation) too difficult to solve. In this 
case other approaches and different estimators have to be chosen. 

An important role is played by the likelihood of Gaussian random variables. If a given observed variable 0„ is 
characterized by gaussianly distributed errors, then it is easy to see that its likelihood is 

£ = e^''\ (A.7) 
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where the chf statistic is defined as: 



X 



(A.8) 



where O'^'"' are the measured values of our observable. In the previous equation we made Oa dependent on a vector 
of parameters A, that we want to fit. Our observable cold be for example the CMB angular power spectrum Cc, 
the primordial power spectrum P(k) or, like in our case the angular bispectrum Bc^^c^j^, and we might be interested 
in knowing the sensitivity of our observation to any cosmological parameter. Our statistical estimate of A will be 
obtained by minimizing That is clearly equivalent to maximize the likelihood. Let us now for simplicity work in 
the one-dimensional case (i.e. our observable depends on a single parameter) and expand the;^^^ about its minimum, 
that is about the best fit value of the parameter A 



X\A)^X\A) + 



1 



2 



(A - AY 



(A.9) 



The linear term vanishes here since we are in the minimum. The quadratic term represents the curvature and defines 
the error on A. If the;^'^ moves very quickly away from its minimum, then our determination of A will be more precise, 
while the error on A will be much larger otherwise. If we define 



1 d'^x^ 
F = — 

2 dA^ 



(A. 10) 



then we can estimate the minimum possible error on /I as 1 / Vf'- It is easy to see that the curvature of the likelihood in 
the Gaussian case matches exactly the definition of Fisher matrix given above. The 1 fsqrtF lower limit on the error 
bar then coincides, as it should, with the Rao-Cramer bound. This at the same time validates the choice of 1 / Vf as 
the error on the parameter, and also shows a simple way to interpret the Rao-Cramer bound. Since the Fisher matrix 
represents the curvature of the In of the likelihood around its maximum, it also provide an intrinsic minimum error 
on the measurement of the parameter A likelihood strongly peaked around its maximum for a given parameter will 
provide stronger constraints on that parameter and vice- vers a. We have however to keep in mind that the curvature F 
constructed above is the curvature of the likelihood only if the distribution of our observable Oa is Gaussian. This, 
strictly speaking, is in general not true, but it is a reasonably good approximation in most cases^". The Fisher matrix 
for any observable is then defined as the second derivative of the;^^^ statistic (|A.8|l. If we compute it explicitly we get 



= Z 



i^OaY 



dOa 

dA 



d^Oa 



dA^ 



(A.ll) 



The se cond term in the sum above is generally neglected. The idea, as explained in |Dodelson| ( |2003] l or in |Press et al.\ 
( I992I is that the observed Oa will oscillate around their real value, making the difiference [Oa - O"*') oscillate around 



zero, resulting in cancellations. We are then left with the expression generally used in the literature, 

21 



1 



i^OaY 



dOa 

dA 



(A. 12) 



In the review we have applied the basic concepts described in this Appendix to the estimation of the non-Gaussian 
parameter f^^"^^^ from the bispectrum of CMB and LSS datasets. We would like to stress again that we have just 
very quickly sketched some essential concepts in estimation theory here. For excellent and much more comprehensive 
reviews of ideas and applications of estimation theory to cosmology we refer the reader to Dodelson ( 2003| l; Martinez 
et al. (2009 i;|Tegmark et al. ( 1997| l. The brief review provided here was actually largely inspired by these works. A 
detailed and complete book about statistical methods and estimation theory is e.g. ,Kendall & Stuart ( 1979 1. 



A clarifying example is provided by the CMB angular power spectrum. We know that C[ is distributed like a^^ with 2/ + 1 degrees of freedom, 
that rapidly gets close to a Gaussian as t grows. 
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